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Abstract 



Weakly Interacting Massive Particles (WIMPs) are one of the leading candidates for 
Dark Matter. Currently, the most promising method to detect many different WIMP 
candidates is the direct detection of the recoil energy deposited in a low-background 
laboratory detector due to elastic WIMP-nucleus scattering. So far the usual procedure 
has been to predict the event rate of direct detection of WIMPs based on some model(s) 
of the Galactic halo from cosmology and of WIMPs from elementary particle physics. 

The aim of this work is to invert this process. In this thesis I will present methods 
which allow to extract information on the WIMP velocity distribution function as well 
as on the WIMP mass from the recoil energy spectrum as well as from experimental data 
directly. 

At first I will derive the expression that allow to reconstruct the normalized one- 
dimensional velocity distribution function of WIMPs from the recoil spectrum. I will also 
derive the formulae for determining the moments of the velocity distribution function. All 
these expressions are independent of the as yet unknown WIMP density near the Earth 
as well as of the WIMP-nucleus cross section. The only information about the nature of 
WIMPs which one needs is the WIMP mass. 

Then I will present methods that allow to apply the expressions directly to exper- 
imental data, without the need to fit the recoil spectrum to a functional form. These 
methods are independent of the Galactic halo model. The reconstruction of the velocity 
distribution function will be further extended to take into account the annual modulation 
of the event rate. 

Moreover, I will present a method for reconstructing the amplitude of the annual mod- 
ulation of the velocity distribution. The only information which one needs is the measured 
recoil energies and their measuring times. An alternative, better way for confirming the 
annual modulation of the event rate will also be given. 

Finally, I will present a method for determining the WIMP mass by combining two (or 
more) experiments with different detector materials. This method is not only independent 
of the model of Galactic halo but also of that of WIMPs. In addition, some meaningful 
information on the WIMP mass can already be extracted from less than one hundred 
events. 
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Chapter 1 
Dark Matter 



One of the most fundamental open questions in cosmology and elementary particle 
physics today is what is the nature of Dark Matter. Earlier the question was whether 
Dark Matter actually exists. But nowadays we have some strong evidence to believe that 
something which we do not know exists. 

As introduction I review briefly the history of the discovery of (the existence of) Dark 
Matter in the Universe. It will be seen that, according to some astronomical observations 
and measurements, more than 80% of the total mass content of the Universe consists of 
Dark Matter. I will also present some models of Dark Matter halo in this chapter. 

1.1 Evidence for Dark Matter 

We call such something "dark" because it (almost) neither emits nor absorbs elec- 
tromagnetic radiation. Historically the observational evidence for the existence of Dark 
Matter came only from galactic dynamics and are gravitational [1] . The following discus- 
sions in this section show that the observed luminous objects (stars, gas clouds, globular 
clusters, or even entire galaxies) can not have enough mass to support the observed grav- 
itational effects [1]. 

1.1.1 Clusters of galcixies 

Clusters of galaxies are the largest gravitationally-bound objects in the Universe. For 
example, our Milky Way and the M31 galaxy belong to the Local Group of Galaxies and 
are part of the Virgo Superclustcr of Galaxies. 

At the beginning of the 1930s, F. Zwicky and other astronomers measured the total 
mass of a few clusters of galaxies and the masses of the luminous objects in these clusters 
of galaxies [2], [3]. Their measurements showed that the masses of these clusters of 
galaxies required to gravitationally bind their galaxies are much larger than the sum of 
the luminous masses of their individual galaxies. 
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Figure 1.1: Rotation curve for OB associations in M31, as a function of distance from 
the galaxy center (figure from 



1.1.2 Rotation curves of spiral galaxies 

The most convincing evidence for the existence of Dark Matter came from the mea- 
surement of the rotation curves of spiral galaxies in the 1970s by V. C. Rubin and other 
astronomers [4]- [6]. 

According to Newton's Second Law, the rotational velocity v of an object on a stable 
orbit with radius r from the center of galaxy is ^ 

v\r) GNM{r) 



namely, 



vlr] oc 



/M(r) 



(1.1) 



:i.2) 



where M(r) is the mass inside the orbit. For an object outside the visible part of the 
galaxy, one would expect that 

1 



v(r] oc 



:i.3) 



^Here the galaxy is assumed to be spherical symmetric. 
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Figure 1.2a: Rotation curves for some simple spiral galaxies. The rotation curves of the 
individual components: visible component, gas, and dark halo, are also shown (figure 
from [7]). 
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Figure 1.2b: Rotation curves for some lower luminosity galaxies. The rotation curves of 
the individual components: visible component, gas, and dark halo, are also shown (figure 
from [7]). 



However, measurements of the circular velocities of clouds of neutral hydrogen in galaxies 
by using their 21-cm emission [1] showed that the rotation curves of spiral galaxies are 
flat (see Figs. 1.1 and 1.2a) or even rising (sec Fig. 1.2b) at distances far away from their 
stellar and gaseous components [4]-[8]. This implies the existence of a "dark halo" around 
the galaxy with a total mass profile: 

M(r)(xr, (1.4) 

i.e., the profile of the mass density should be 

p(r)oc^, (1.5) 
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Figure 1.3: The position of the Sun in the Milky Way. The visible (luminous) component 
has been shown. It can be seen that our Solar system is already out of the Central Bulge 
of the Galaxy. 



smce 

rr 

„'2 



M(r) = 47r / r'^p{r')dr' . (1.6) 







1.1.3 Escape velocity from the Milky Way 

The escape velocity from the Milky Way at the position of our Solar system has been 
estimated as [9], [1] 

^et''^"> 450 km/s. (1.7) 

It is much larger than can he accounted for by the luminous matter in our Galaxy. It 
is not difficult to understand why this result so surprising if one thinks about the huge 
difference between the escape velocity from the Sun's surface [1]: 

= 617.5 km/s , (1.8) 

and that from the Solar system at the position of the Earth [10]: 

C^ = 42.1km/s. (1.9) 

Recall that the gravitational well in our Solar system is essentially only caused by the 
Sun's mass which dominates the total mass of the Solar system. If the mass of the 
luminous matter in our Galaxy would also dominate the total mass of the Galaxy (see 
Fig. 1.3), the escape velocity from our Galaxy at the position of our Solar system would 
also be reduced (at least) one order of magnitude. 
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1.2 Cosmological density parameters 

The cosmological density parameter of a given component of the total energy of the 
Universe i has been defined as the density of this component averaged over the Universe, 
Pi, in units of the critical energy density of the Universe, pcrit, 

Qi^-^. (1.10) 

Pcrit 

The critical energy density of the Universe is the value that makes the geometry of the 
Universe fiat (a more detailed explanation about the "flat Universe" will be given in 
Subsec. 1.2.2) [11]: ^ 

Pcrit = 4^ - 2.775 X Mq/Mpc' 

~ 1.878 X IQ-^^h^ g/cia . (1.11) 

Here H is the expansion rate of the Universe (the time dependent Hubble parameter) 
defined as 

H=^ (1.12) 
d 

with the scale factor of the Universe, a{t), and Hq denotes the expansion rate of the 
Universe "at the present epoch" (redshift z — 0), 

Ho = lOO/i km/s/Mpc , (1.13) 

with the dimensionless Hubble constant h. Moreover, the Newtonian gravitational con- 
stant, the mass of the Sun, and the parsec (pc) are given as [11] 

Gn = 6.674 X 10-^^ mVkg/s^ , (1.14) 

Mo = 1.988 X 10^° kg, (1.15) 
1 AU 

1 pc = ~ 3.0857 X 10^*^ m ~ 3.262 c ■ yr , (1.16) 

1 arc sec 

where the astronomical unit (AU), i.e., the mean distance between the Earth and the 
Sun, and the speed of fight, c, are given as [11] 

1 AU = 1.4960 X 10^^ m, (1.17) 

c = 2.99792458 x 10^ m/s . (1.18) 

In the rest of this section I present briefiy some important astronomical measurements 
and their current results, by which the cosmological density parameters of different com- 
ponents of our Universe can be determined pretty exactly (to one or even two significant 



^Note that pcrit here is the critical "energy density". However, the factor in the expression has 
been usually neglected. 
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figure accuracy [12]). Particularly prominent are the measurement of the anisotropy of 
the cosmic microwave background (CMB) radiation, led by the three-year results from 
the Wilkinson Microwave Anisotropy Probe (WMAP) [13]. In the last subsection we will 
see that the cosmological density parameters also show the evidence for (or the necessity 
of) the existence of Dark Matter (and, more exactly, also of Dark Energy, both of them 
are "something which we do not know"). More details about theoretical explanations 
and experimental results of these measurements can be found in e.g., Refs. [14], [15], [12], 
and [16]- [18]. 

1.2.1 Cosmic microwave background (CMB) 

The cosmic microwave background radiation (CMBR or CBR) discovered in 1965 
provides the fundamental evidence for the hot Big-Bang model of the early Universe [14]. 
The spectrum of the CBR can be described very well by a blackbody function with the 
temperature T [18]. The energy density of "CMB photons" can then be obtained directly 
as [11] 

p = (1 19) 

15 {her ' ^ ^ ^ 

where the Boltzmann's constant, ks, and the Planck's constant, h, have been given as 
[14] 

kB = 1.3807 X 10-23 J/K , (1.20) 

^1 = 1.0546 X 10-=^^ Js. (1.21) 

The present (mean) CBR temperature has been measured as [11] 

To = (2.725 ±0.001) K. (1.22) 

1.2.2 Anisotropy of the CMB radiation 

Another important observable quantity from the CMB is its anisotropy: tiny tem- 
perature difference (of order of 10"^ of the magnitude of the mean temperature T [18]) 
between two points on the sky (see Fig. 1.4). The measurement of the anisotropy of the 
CBR can be expanded in spherical harmonics as: 

5T{e,<f>)^J2^lmyim{d,<P). (1.23) 
l,m 

Here the multipole number I is given as 
200° 

/ - ^ , (1.24) 
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Figure 1.4: Anisotropy of the CMB radiation. The detailed, all-sky picture of the in- 
fant Universe from three years of WMAP data. The image reveals 13.7 billion year old 
temperature fluctuations (shown as color differences) which correspond to the seeds that 
grew to become the galaxies. The signal from our Galaxy was subtracted using the 
multi-frequency data. This image shows a temperature range of ±200 fiK (figure from 
NASA/WMAP Science Team). 



and a useful quantity Ci has been defined as 

Q^{hm?) = -^ E \ai^\\ (1.25) 

m=—l 

The anisotropy of the GBR offers the best means for determining the curvature of the 
Universe, -Rcurv? [15] and thereby the "total matter/energy density" of the Universe, Qq, 
according to the Friedmann equation [15], [12]: 

^0-1 = ^5^, (1.26) 

where is a curvature constant which can be chosen to take only three discrete values: ±1 
and 0. According to the Friedmann equation, when the total matter/energy density of our 
Universe is equal to 1, the Universe is "spatially flat" (-Rcurv = oo, or, equivalently, k = 0). 
While, for > 1 (^o < 1); the constant k should be +1 (—1) and we call the Universe 
"closed" ("open") [16]. In Fig. 1.5 one can find that the anisotropy power, sometimes 
shown as /(/ + 1) Ci/2tt, oscillates (the so-called "gravity-driven acoustic oscillations") 
with some "acoustic peaks" . Roughly speaking, the angular position of these peaks is a 
sensitive probe of the spatial curvature of the Universe: if our Universe is open (close), 
these peaks should lie at higher (lower) / [18]. 

Moreover, according to standard Big-Bang Gosmology, the higher the primordial mat- 
ter density, the shorter the duration of the epoch of structure formation and thereby the 
larger fluctuations in the GBR [1] , or, equivalently, the stronger these acoustic oscillations 
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Figure 1.5: The angular power spectrum of the CMB temperature from three-year data 
of the WMAP satelhte. The sohd curve is the prediction from the best-fitting ACDM 
model. The error bars on the data points (which are tiny for most of them) indicate 
the observational errors, while the shaded region indicates the statistical uncertainty 
from being able to observe only one microwave sky, known as cosmic variance, which is 
the dominant uncertainty on large angular scales [12]. The first peak around / ~ 200 
corresponds to ~ 1° (figure from NASA/WMAP Science Team). 

[15]. Hence, the relative height of the first acoustic peak can be used to determine the 
"primordial matter density" . 

More details about the physics and the analyses of anisotropy of CBR can be found 
in Ref. [18]. 

1.2.3 Age of the Universe 

As mentioned in the previous subsection, the higher the primordial matter density, the 
faster the Universe expanded and thus the shorter the age of the Universe reaching its 
present size. Hence, the measurements of the age of the Universe, Tu, and the expansion 
rate of the Universe, h, can give the upper and lower limits on the "matter density" in 
the Universe. 

According to WMAP results combined with other astronomical measurements, the 
age of the Universe has been estimated as [13], [11] 

Tu = 13.7l°iGyr. (1.27) 



9 



1.2.4 Present expansion rate of the Universe 

According to the Hubble law [15]: 

/f„=^. (1.28) 

Here the velocity v can be determined by the redshift, thus the most accurate direct 
methods for measuring distances to distant objects d can be used to estimate the Hubble 
parameter Hq [19]. Currently, there are two methods for measuring extra galactic dis- 
tances [19]: time delays between luminosity variations in different gravitationally lensed 
images of distant quasars and the Sunyaev-Zel'dovich effect: Compton scattering of the 
CMB by the hot electrons in clusters of galaxies. Note that the error on the estimates 
of the Hubble parameter is dominated by one systematic uncertainty: the distance from 
out Galaxy to the Large Magellanic Cloud (LMC), which has been used to cahbrate the 
Cepheid period-luminosity relationship [19]. 

The dimensionless Hubble constant has been estimated as [13], [11] 

/i = 0.73ini, (1-29) 
and the present expansion rate of the Universe can then be given as 

Hq = 73tl km/s/Mpc . (1.30) 



1.2.5 Abundances of the light elements 

BBN predicts the primordial abundances of the light elements. Thus measurements 
of the primordial abundances of the light elements produced in the Big Bang, such as 
deuterium (D), helium (^He and "^He), and lithium (^Li), can also give the upper and 
lower limits of the "baryon density" in the Universe. 

Moreover, among these four light elements, because the primordial abundance of 
deuterium depends strongly on the baryon density (oc p^^''^), and it can only be destroyed 
by the astrophysical processes, deuterium becomes the most powerful "baryometer" [15]. 

Fig. 1.6 shows the theoretically predicted abundances of the four hghtest elements 
and the observational results. 



1.2.6 Gas- to- total mass ratio 

The clusters of galaxies formed due to density perturbations with a co-moving size 
of the order of 10 Mpc and gathered material from such a large region of the Universe 
[15]. Meanwhile, most of the baryons in the clusters of galaxies reside not in the galaxies 
themselves but in form of hot intercluster, x-ray emitting gas [15]. Hence, by measuring 
the gas-to-total mass ratio of the cluster, /gas/totab and combining with the (measured) 
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Figure 1.6: The predicted abundances of ^He (mass fraction), D, ^He, and ^Li (number 
relative to hydrogen) by the standard model of the BBN as a function of the baryon 
density. Yp = 2n^n-p/{ny^ + rip) ~ 0.25, where and n-p are the neutron and proton 
number densities. Widths of the curves indicate 2a theoretical uncertainty. Boxes indi- 
cate the observed light element abundances (smaller boxes: 2cr statistical errors; larger 
boxes: ±2o" statistical and systematic errors). The narrow vertical band indicates the 
CMB measurement of the cosmic baryon density (figure from [17]). 
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baryon density in the Universe, Q^, we can determine the "matter density" in the Universe 
[15]: 

^lm = -r^. (1.31) 

/gas/total 

There arc two methods for determining the mass of the intcrchistcr gas: the x-ray 
flux emitted from the intcrcluster gas or the Sunyaev-Zel'dovich CBR distortion caused 
by CMB photons scattering off hot electrons in the intracluster gas [15]. While, there 
are also three independent methods for estimating the total mass of a cluster: the mo- 
tions of cluster galaxies with the virial theorem, assuming that the gas is in hydrostatic 
equilibrium and using it to infer the underlying mass distribution, or mapping the cluster 
mass directly by gravitational lensing [15], [20]. Within their uncertainties and where 
comparisons can be made, the two methods for determining the mass of the intercluster 
gas and the three methods for estimating the total mass of a cluster are consistent with 
each other, respectively [15]. 

1.2.7 Mass-to-light ratio 

One other method for estimating the "matter density" of the Universe is using the 
mass-to-hght ratios [15]: 

Pm^[j^)C, (1.32) 

where C is the averaged luminosity density of the Universe [1], [15]. In V-band [1] and 
in B-band [15], we have, respectively, 

Cy = (1.7 ± 0.6) X lO^h Lo/Mpc^ , (1.33a) 

and 

jCb = 2.4 X lO^h Lo/Mpc^ , (1.33b) 

where Lq is the luminosity of the Sun [11], 

Lq = (3.846 ± 0.008) x 10^^ W . (1.34) 

Once we have estimated the mass-to-light ratios of some systems, i.e., 
M 

T, = — , x = V, B. (1.35) 

Then, combining Eqs.(l.lO), (1.11), (1-32) and (1.33a) or (1.33b), the total matter density 
can be obtained as 
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Here 



Cv = 6.1 , 



Cb = 8.6, 



(1.37) 



and T, 



o is the mass-to-light ratio of the Sun, 



To = 5.169 X 10^ kg/W, 



(1.38) 



where I have used Eqs.(1.15) and (1.34). 

1.2.8 Supernovae type la (SNe la) at high-redshift 

If we could measure the present extra galactic distances do and velocities vq, they 
should obey the Hubble law [15]: 



since the expansion of the Universe is simply a rcscaling. But what we can actually 
measure are the distances dz and velocities Vz at an earlier time (redshift z). If we 
suppose that the expansion of our Universe should slow down due to the attractive force 
of gravity, i.e., Hz > Hq, the measured galactic velocities should be larger than that 
expected by the Hubble law: 



or, equivalently, for the galaxies with known velocities, their distances should be shorter 
than that expected by the Hubble law: 



In 1998 two groups: the Supernova Cosmology Project and the High-z Supernova 
Search Team have pubhshed their "magnitude-redshift (Hubble) diagram for fifty-some 
type la supernovae (SNe la) out to redshifts of nearly 1" [15]. By using SNe la as standard 
candles for estimating the distances to faraway galaxies, the two groups concluded that 
distant galaxies are moving slower than predicted by the Hubble law and that this implies 
an accelerated expansion of our Universe [15]. 

In order to explain this observational indication, i.e., in order to find the discrepancy 
between the measured total matter (/energy) density, Qq, and the matter density, Qm, 
(data given in the next subsection), a new term "Dark Energy" ^ (or sometimes also 
called 'Viuintessence" ) has l^een introduced [15]. 

^Dark Energy is beyond the area of my research and thus will not be discussed in this work. Short 
reviews and summaries can be found in e.g., [15], [12], [16], and [21]. 




(1.39) 



Vz = dzHz > dzHo , 



(1.40) 
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1.2.9 Cosmological densities of different components 

According to the various astronomical measurements described above (and other 
measurements, e.g., the pecuhar velocities of galaxies, the shape of the present power 
spectrum of density perturbations, and the opacity of the Lyman-a forest toward high- 
redshift quasars), we can conclude today the cosmological densities of different compo- 
nents as follows. 

The total matter/energy density is [13], [11] 

Oo = 1.003l°:°i?. (1-42) 

It can be separated into Dark Energy [11]: 

^^A^^ = 0.76^^6, (1-43) 
where we have used 

and total matter [13], [11]: 

n^ = 0.127t'oZh-' = Q.24.t'o'ol (1-45) 
The total matter in the Universe can also be separated into baryons [13], [11]: 

Q, = Om23troZh-' = OM2t',Z,, (1.46) 
and non-baryonic Dark Matter [11]: 

^^DM = Q„ - = 0.105iH?o = 0.20iHI ■ (1-47) 
Among the baryons in the Universe there is luminous matter with a density of [1] [22]: 

aum^O.Ol, (1.48) 
including the density of the stars [23] : 

^^stars = (0.0023 ~ 0.0041) ± 0.0004 . (1.49) 

While, the non-baryonic Dark Matter can be separated into Cold Dark Matter (CDM) 
and Hot Dark Matter (HDM) (the definitions and some discussions about the CDM and 
HDM will be given in Sees. 2.1 and 2.2, respectively). Finally, among the relativistic 
particles (see Sec. 2.2), the density of the CMB photons can be estimated directly by 
inserting Tq in Eq.(1.22) into Eq.(1.19) as [11] 

= (2.471 ± 0.004) X 10"^ h-^ = (4.6 ± 0.5) x 10"^ , (1.50) 

and the density of the neutrinos has been estimated as [11] 

< 0.014 (95% C. L.) . (1.51) 
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1.3 Galactic halo models 



In this section I present some simple halo models. For estimating some characteristics 
of halo Dark Matter, such as the velocity dispersion of Dark Matter particles, -u, and the 
local Dark Matter density, po, the rotation curve of our Galaxy is the most important 
observational quantity, since it measures the change in density and sets the scale for the 
depth of the Galactic potential well [1]. Essentially, all important direct information 
which has been obtained about the halo is provided by the rotation curve [24]. However, 
due to our location inside the Milky Way, it is more difficult to measure the accurate 
rotation curve of our own Galaxy than those of other galaxies (see Fig. 1.7). In addition, 
the "disk contribution" to the rotation curve must be known to infer the halo contribution, 
but precise determination of the disk contribution is also difficult. 

1.3.1 Standard assumptions of Dark Matter halo 

The velocity dispersion of Dark Matter particles ^ in the Solar neighborhood has been 
assumed as 

= (t;2)V2 ^ 270 km/s. (1.52) 

And the lAU standard value for the rotational velocity at the Sun's distance from the 
Galactic center is [11] 

vo = v,,t{ro) ~ (220 ± 20) km/s , (1.53) 

where the distance from the Sun to the Galactic center is [11] 

ro ~ (8.0±0.5) kpc. (1.54) 

On the other hand, the local Dark Matter density (the Dark Matter density near the 
Solar system) is given by [1] 

Po = p{ro) f« 0.3 GeV/cVcm^ 

5 X 10"^^ g/cm^ (1.55) 

with an uncertainty of slightly less than a factor of 2 [1], [15]. Here I have used [25] 

1 GeV/c^ = 1.7827 x 10"^^ g . (1.56) 



^ Strictly speaking, the velocity dispersion should be (v^) — (v)^. However, since the major component 

of Dark Matter should be cold (a detailed discussion will be given in Sec. 2.1) and assumed to have 
negligibly small velocity average (v) in the Galactic rest frame, its rms velocity {v'^) has been called 
sometimes simply as the velocity dispersion. 
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(]a]£ctocenirLc Radius | kpc | 



Fi^rfr 2, Thero^anion curve for 4he Milky Wiiy for vjIik of — 7. 1 kpc. 

= Ifi5km£-', and = fi.^kpc. =^a}l;m£-'. The figure ilso 
diQW£ one cf ilK in which the rolJilicin curve can bedcccmpDscd i nio 
the com ribm ions Tie 111 dirreren* mass coiifionenls: 4hebul^ idolted line): 
the sidlardisc inilcd circl cs): 4 he H I layer IcrosseE, where ne^ajl i vc values 
mean than 4he fcrce is dimmed outwaids): 4he layer iclicles): and 4he 
dark halo ida^ed line]. The best-dMln^ model, which is cbtaincd by 
summing 4he individual CDmpcnents i n quadrdire, is shown as a full line. 



Figure 1.7: Rotation curve for the Milky Way, as a function of distance from the Galactic 
center, with two different assumptions for the Sun's distance from the Galactic center. 
To, and the rotation velocity at ro, Vroti^o) (figure from [8]). 
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1.3.2 Canonical isothermal spherical halo model 



The simplest halo model is an isothermal spherical halo. An empirically plausible 
radial profile for a spherical galactic halo is constrained only by its contribution to the 
galactic rotation curve. This means that the radial profile should approach to a constant 
near its core so that it gives rise to a linearly rising rotation curve at small radii, and it 
should fall as or eventually faster at large radii to provide a fiat rotation curve [24] . 
The density profile of the cored isothermal spherical halo is given by [1] 



Pis{r) = Po 



(1.57) 



where po is the local halo density and is the core radius of the isothermal spherical 
halo, within which the density Pis{t) behaves no longer as 1/r^, but goes to a constant 
as r approaches 0. 

Substituting this expression into Eq.(1.6) and using Eq.(l.l), the rotational velocity 
at a radius r from the halo center can be found as 



V) = AttGn ■ - / r'>(r') dr' 



r Jo 



1 _ - tan-^ - 



(1.58) 



where I have used 
dx 



I 



X 



a tan 



+ x'^ 

Define as the (measured) rotational velocity as r — > oo. One can find that 

= Visir ^ oo) = AttGnPo {r^ + ^o) > 
thus the local halo density po in Eq.(1.57) can be expressed as 



Po 



;i.59) 



;i.60) 



Meanwhile, combining Eqs.(1.58) and (1.59), the core radius of the isothermal spherical 
halo in unit of ro, i.e., Tc/ro, can be solved (numerically) by [1]: 



— I tan 



^1=1- 



(1.61) 



Eqs.(1.60) and (1.61) show how we can estimate the local halo density and the halo core 
radius in the isothermal spherical halo model once the rotational velocities Vq and 
have been measured. 

Finally, substituting Eq.(1.60) into Eqs.(1.57) and (1.58) the density profile and the 
rotation curve of the isothermal spherical halo model can be rewritten as 

o,2 



Pis(r) 



1 



(1.57') 
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Figure 1.8: (a) The radial density profile pis{r) given in Eq.(1.57') and (b) the rotation 
curve fis('") given in Eq.(1.58') of the canonical isothermal spherical halo model. Here I 
have used Voo = 220 km/s [24] and = 2.6 kpc (see Subsec. 1.3.4). 
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and 



vis{r) 



1/2 



:i.58') 



Figs. 1.8 show the radial density profile (upper frame) and the rotation curve (lower 
frame) of the canonical isothermal spherical halo model given in Eqs.(1.57') and (1.58'), 
respectively. 

1.3.3 Alternative isothermal spherical halo model 

An alternative density profile for the isothermal spherical halo has been given by [24] 



PAis{r) = po 



;i.62) 



Tc + r 

where po is again the local halo density and Vc is the core radius of this alternative 
isothermal spherical halo mode. ^ Using Eqs.(1.6) and (1.1), the rotation curve of this 
alternative isothermal spherical halo can be found as 



^Ais('^) = 47rGivPo {tc + roY 



1 + 



where I have used 
dx 



I 



1 



{ax + by 



{ax + b) — 



In 



- 2b In (ax + b) 



(1.63) 



ax + b 

Meanwhile, the core radius of this alternative halo model in unit of Tq can be solved 
(numerically) by : 



2ar In 



«c + 1 



1 - 



^Ais('^o) 



(1.64) 
(1.65) 



etc / etc + 1 

where I have defined 

_ fc 
ac = — . 

Finally, from Eq.(1.63), one has 

= ^Ais(^ ^ oo) = A-kGnPq {tc + rof , (1.66) 

the density profile and the rotation curve of the alternative isothermal spherical halo 
model in Eqs.(1.62) and (1.63) can be rewritten as 



47rGjv Vtc + r 



and 



^^Ais('^) = 



1 + 



- 2 



In 



rc + r 



1/2 



fl.62' 



(1.63') 



Tc + r \ r y \ Tc 
Fig. 1.9 shows the equations for solving the ratios of the core radii of two cored 
isothermal spherical halo models to the distance from the Sun to the Galactic center, 
ttc = Tc/tq, given in Eqs.(1.61) and (1.64). 



'"'Here, for simplicity, I use the same notation as in Eq.(1.57), but Vc for these two halo models are 
not the same. The determination of these core radii will be given in Subsec. 1.3.4 
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alpha_c 



o.s 



Figure 1.9: The equations for solving the ratios of the core radii of two cored isothermal 
spherical halo models to the distance from the Sun to the Galactic center, ac = Tc/to- 
The solid (red) line shows Eq.(1.61) and the dash-dotted (blue) line shows Eq.(1.64). 
Here the dash (green) line denotes fis(?"o) = VAisi^o) = fhaio('"o) = 170 km/s [1], [24] and 
Voo = 220 km/s. 

1.3.4 Evans' power- law halo model 

Even if we consider only spherical halo distributions, there are still some latitudes in 
our choice for the precise form of the radial density profile of a halo model. Meanwhile, 
N-body simulations of gravitational collapse produce axisymmetric or triaxial halos [26], 
and other spiral galaxies appear to have flattened halos [27], [26]. 

The equipotentials of elliptical galaxies and the halos of spiral galaxies could be 
roughly stratified on similar concentric spheroids [28]. This suggests a useful approx- 
imation to their gravity field as [28]: 

where ipa is the central potential, Tc is the core radius of the halo, and q is the axis 
ratio of the equipotentials or the so-called flattening parameter. Note that I use here the 
cylindrical coordinates and r denotes the distance from the point which one considers 
to the rotation axis. The potential is just a power of the spheroidal radius r, thus this 
model has been called the power-law halo model [28]. 
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Using the Poisson Equation: ^ 

= -4:TiGnP, 



;i.68) 



the density distribution of the gravitational potential described in Eq.(1.67) can be ob- 
tained as [28] 



p{r, z) 



{2q^ + l)rl + (1 - (3qy^ + {2q^ -(3- l)z'/q^ 
{rl + r-^ + z^/q^Y/'^+-^ 



where 



:i.69) 



(1.70) 



Meanwhile, the velocity of the circular orbit in the equatorial plane with radius r can be 
obtained as [28] 



^^circ('^) = Va 



c 



1/2 



(^2 _^ ^2j/3/2+l 

since the central force 

Fcen{r, z = 0) = -Vipir, z) ■ r 



^circ('^) 



z=0 



(1.71) 



;i.72) 



The rotational velocity in Eq.(1.71) is asymptotically falling if /3 > and rising if /3 < 
[28]. 

On the other hand, the model with spheroidal equipotentials and a completely flat 
rotation curve at large radii is well known as an axisymmetric logarithmic potential [29], 
[30], [28]: 



y2 / ^2\ 



;i.73) 



where Va is the rotational velocity at large radii (i.e., Voo used in the previous two sub- 
sections) . Using the Poisson Equation, the density distribution can be found as [30] 



p(r, z) 



(1.74) 



47rGAr [ q'^{r'^ + r'^ + z'^/q'^y 

Comparing this expression with the expression in Eq.(1.69), it can be seen that the 
logarithmic potential given in Eq.(1.73) has the properties corresponding to the missing 
(3 = case in Eq.(1.67) [28]. The velocity of the circular orbit in the equatorial plane 
with radius r due to the potential given in Eq.(1.73) can be obtained as 

/ ^2 \ 1/2 



^^circ('^) = Va 



J.I _|_ 



(1.75) 



^Note that the Poisson Equation holds actually for the "total" potential and the "total" density 
distribution, i.e., for luminous baryonic matter and Dark Matter together, not for each component 

separately. 
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Furthermore, let q — 1 and use the replacement: 



2,2 2 

r + z ^ r , 



;i.76) 



we can rewrite the potential in Eq.(1.73) to the potential for a spherical halo as [30] 



^(r) = -^In (r^ + r^) 



;i.77) 



The density distribution and the velocity of the circular orbit of this spherical Evans 
model with radius r can be obtained as, respectively [30], [24], 



PPL[r) 



AttGn 



{r1 + r'^Y 



and 



1/2 



rf. + 



;i.78) 



(1.79) 



As done in Subsecs. 1.3.2 and 1.3.3, the core radius of the halo in unit of Tq can be solved 
by means of the following equation [24] : 

1/2 

- 1 . (1.80) 



Finally, using the local rotational velocity given in Eq.(1.53): 

^;,ot(ro) :^ 220 km/s, (1.53') 
and the assumption for the disk contribution to the rotational velocity at r = Tq [1], [24]: 

^^disk(ro) = 140km/s, (1.81) 
the local halo contribution can be found as 

t'haio(ro) ^ 170km/s, (1.82) 

since, 

^rot(ro) = \/^^disk('^o) +^^halo('^o)- (1-83) 

Then the core radii for the canonical, alternative, and Evans spherical halo models given 
in Eqs.(1.57'), (1-62'), and (1.78) can be found by means of Eqs.(1.61), (1.64), and (1.80) 

as 



rc,is 2.6 kpc , 
rc,Ais - 0.86 kpc , 



(1.84) 

(1.85) 



and 



^cPL - 6.6 kpc , 



;i.86) 
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r[kpc] 

(b) 

Figure 1.10: (a) The radial density profiles and (b) the rotation curves of different isother- 
mal spherical halo models. The solid (red) lines indicate the canonical halo model in 
Eq.(1.57') with Tc ~ 2.6 kpc, the dash (blue) lines indicate the alternative halo model in 
Eq.(1.62') with Vc — 0.86 kpc, and the dash-dotted (black) lines indicate the spherical 
Evans' halo model in Eq.(1.78) with rc — 6.6 kpc. Here I have used foo = 220 km/s [24]. 
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where I used 



ro~8.0kpc, (1.54') 

and the rotational velocity at large radii as [24] 

v^ = Va = 22Q\im/s. (1.87) 

Figs. 1.10 show the radial density profiles (upper frame) and the rotation curves (lower 
frame) of different isothermal spherical halo models with the core radii obtained above. 

1.3.5 NFW halo model 

Besides the three spherical halo models presented in the previous three subsections, 
there is also the well-known NFW density profile given as [31] -[33]: 

Pnfw(?^) = 1—, — rjT^. — -, — ■ (1- 

Here ps is the characteristic density [33] : 



(1 89) 

_ln(l + c,,)-c,/(l + c,)J ' ^ ■ ^ 

with the virial overdensity and the halo concentration parameter Cg, and is the 
characteristic radius [33]: 

r^ir 1.2 X 10^ / Mv 



' vir 




The virial radius, Tyir, defined as the radius, inside which the average overdensity is 
times the critical density of the Universe, and the virial mass, Myi^ is then the total mass 
within this virial radius Tyir [32]. Usually the average overdensity A^ has been chosen as 
200, and the virial radius and the virial mass for As = 200 have been then specially 
labeled as r2oo and M200, respectively [32]. However, for a fiat Universe, i.e., flo — 1, the 
average overdensity A^ has been found to be [34] 

As = IStt^ ~ 178 . (1.91) 

Note that there is a good correlation between and M^ir in Eq.(1.90), which results from 
the fact that dark halo densities refiect the density of the Universe at the epoch of their 
formation and that halos of a given mass are preferentially assembled over a narrow range 
of redshifts [32] . Hence, as lower mass halos form earlier, at times when the Universe was 
significantly denser, they are more centrally concentrated [32]. 

Although the NFW density profile given in Eq.(1.88) doesn't approach the 1/r^ form 
for large r, high-resolution N-body simulations of structure formation and some observa- 
tional results have shown that the NFW profile indeed provides a good description of the 
density distribution in clusters [34], 
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Moreover, the NFW density profile given in Eq.(1.88) has been expanded to the 
following form [35] : 

PNFWE(r) = . , .J[\ I ■ (1-92) 

For the original NFW profile, one has a = 1, = 2; for a modified NFW profile, one can 
use a — V — ?>I2\ and for a so-called Hernquist profile one has a = 1, v — ?> [36], [35]. 
On the other hand, A. Burkert has suggested a Burkert density profile [32]: 

=(1 + [l + (r/n)r 

Some analyses show that density profiles of dwarfs and low surface brightness galaxies 
can be fitted by the Burkert profile much better than the NFW profile [34] . 

1.3.6 Bulk rotation 

So far the halo models presented above arc only classified by their density profiles or the 
gravity field, namely, their mass distribution. Because, as mentioned in the beginning 
of this section, what has been measured is just the rotation curve [24]. However, the 
rotation curve is determined by the halo mass distribution and is insensitive to its velocity 
distribution [24]. Thus, even though there are some theoretical arguments against a 
rotation-dominated velocity distribution, there is no empirical evidence to rule out a 
halo with bulk rotation [24]. Note that such bulk rotation can also affect the velocity 
distribution of WIMPs seen near the Earth [24]. 

Halo models with bulk rotation can be constructed by taking linear combination of 
the velocity distribution function [26] : 

/rot(v) = ai.ot/+(v) + (1 - arot)/-(v) 

arot/(v), for^;^>0, 
= < (1.94) 

(1 - arot)/(v) , for < . 

A non-rotating halo has Orot = 0.5, whereas a counter-rotating or a co-rotating one has 
< arot < 0.5 or 0.5 < Crot < 1, respectively [24], [26]. 

Moreover, Orot is related to a dimensionless spin parameter Arot, which usually has 
been used to quantify the galactic angular momentum [24]: 

A,ot = 0.36|a,ot-0.5| . (1.95) 

Numerical studies of galaxy formation find that |Arot| < 0.05, corresponding to 0.36 < 
a,ot < 0.64 [26]. 



^The velocity distribution of WIMPs will be discussed in Chap. 3. 
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Chapter 2 

Candidates for Dark Matter 



As defined in the Introduction, Dark Matter (almost) neither emits nor absorbs 
electromagnetic radiation. It is thus non-luminous. 

Meanwhile, as described in Sees. 1.1 and 1.2, so far we can "observe" (or, actually, 
"feel") the existence of Dark Matter only through its gravitational effects. Moreover, 
according to the observational results for the rotation curves of spiral galaxies (described 
in Subsec. 1.1.2), Dark Matter forms halos with an approximately spherical distribution 
around galaxies. Hence, Dark Matter (almost) does not interact with ordinary matter and 
is collisionless. Otherwise, if Dark Matter could interact with ordinary matter, it would 
dissipate its kinetic energy after the interactions, fall onto galaxies, settle deep into the 
galactic gravitational wells, and thus form the "galactic disks" with ordinary matter. 

On the other hand, (the major part of) Dark Matter particles should moved non- 
relativistically in the early Universe or, equivalently, have sufficiently low primordial ve- 
locity dispersion, in order to allow it to merge to galactic scale structures (e.g., galaxies 
and clusters of galaxies). In contrast, although neutrinos are also collisionless, they 
moved relativistically and have thus too large velocity dispersion to build galactic scale 
structures. 

Therefore, Dark Matter should be some "non-luminous, non-baryonic, non-relativistic, 
and collisionless" elementary particles have not yet been discovered. In addition, the 
candidates for Dark Matter must satisfy the following cosmological conditions: they 
must be stable on cosmological time scales and have the right relic cosmological density 
[37]. 

In this chapter I present some most motivated and studied candidates for Dark Matter. 
Most of them are non-baryonic and non-relativistic particles. However, some relativistic 
particles and baryonic objects could also be (part of) Dark Matter. 

2.1 Cold Dark Matter (CDM) 

"Cold" has been used here to indicate that such Dark Matter particles moved non- 
relativistically at the matter-radiation decoupling time in the early Universe [22], i.e., at 
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Figure 2.1: Extension of the Standard Model of particle physics to supersymmetric models 

the time in which galaxies could just start to form. Due to their relatively slower veloc- 
ities Cold Dark Matter would first form some relatively small galactic scale structures; 
large galaxies and clusters of galaxies are formed through "hierarchical merging" of these 
smaller structures. 



2.1.1 Minimal Supersymmetric Standard Model (MSSM) 

Supersymmetry has been considered to solve the hierarchy problem in the Standard 
Model (SM) of particle physics: Why is the electroweak scale (-Eew — 0(100 GeV)) so 
small compared to the other known scales such as the grand unification scale (£^gut — 
10^6 GeV [1]) or the Planck scale {Epi ~ 10^^ CeV [1])? 

As shown in Fig. 2.1, supersymmetry provides a natural framework for discussing the- 
ories with large hierarchies of scales and unification with gravity [1]. In supersymmetric 
models, for every fermionic degree of freedom there is a bosonic degree of freedom and vice 
versa. This means also that, for each "normal" particle, there will be a supersymmetric 
partner. Hence, the particle spectrum has been greatly extended in the MSSM. 

The particles of typical supersymmetric models are given in Table 2.1. The spectrum 
of the normal particles is specified in the same manner as in non-supersymmetric models. 
Quark mass matrices determine the masses and the mixing angles, which are encoded 
in the Cabibbo-Kobayashi-Maskawa (CKM) matrix. The pattern of gauge-symmetry 
breaking is unchanged from the Standard Model, and gives the same tree-level relation 
between the masses of the and Z° bosons. Quarks in the Standard Model have spin 
|, while their supcrpartners, squarks, are scalars [1]. There are two squarks (left-hand 
and right-hand) for each quark. In some models there is no mixing between different 
fiavors, and each squark is associated with a given quark [1], for example, ul and ur, dL 
and dji- However, generally the three up-quarks can mix among themselves and similarly 
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Normal particles 


SUSY partners 


Name 


Symbol 


Name 


Symbol 


up-quarks 
down-quarks 


q = u, c, t 
q = d, s, b 


up-squarks 
down-squarks 


UL, UR, CL, CR, tL, f'R 

dL, dR, SL, SR, bL, bR 


leptons 
neutrinos 


e, n, T 


sleptons 
sneutrinos 


CL, eR, flL, fJ'R, TL, TR 


gluons 


9 


gluinos 


9 


photon 
Z boson 


7 


photino 7 
Z-ino Z 


neutralinos 


Xl, X2, Xs, X4 


light scalar Higgs 
heavy scalar Higgs 


hP 
H° 


neutral higgsinos hP , 


pscudoscalar Higgs 


AO 






charged Higgs 
W bosons 




charged higgsinos 
gauginos, W-inos 


charginos 


xf, xt 


graviton 


G 


gravitino 


G 


axion 


a 


axino 


a 



Table 2.1: Particles of typical supersymmetric models 

for the three down-quarks, so there are totally six up-squarks and six down-squarks in 
the particle spectrum [1]. Similarly for the leptons. In these models, left-right sfermion 
mixing is proportional to the corresponding fermion mass [1]. Thus there is little left- 
right mixing for u, d, and s squarks or selectrons or smuons, but mixing of staus and c, 
h, and especially t squarks can be substantial [1]. 

A most important technical difference between the Standard Model and the MSSM 
occurs in the Higgs sector. Two weak isospin Higgs doublet fields arc required in the 
MSSM, whereas only one is required in the SM [1]. This enrichment of the Higgs sector 
gives rise to five physical states and provides an important phenomenological window. 
The superpartners of the and charged Higgs bosons, the gauginos and the charged 
higgsinos, carry the same 5'[/(3) x ^7(1) quantum numbers. Thus they will generally mix 
after electroweak-symmetry breaking, and the two resulting mass eigenstates are linear 
combinations known as charginos [1]. Similarly for the superpartners of the photon, 

boson, and neutral Higgs bosons. These fields generally mix to create four mass 
eigenstates called neutralinos [1]. In many supersymmetric models, constraints on the 
Higgs-higgsino sector are therefore an important area of supersymmetric phenomenology 

[1]- 

In Table 2.1, the tilde ~ has been used to denote a supersymmetric particle. However, 
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the tildes for neutralinos and charginos are sometimes omitted since there is no ambi- 
guity for such particles. Moreover, the lightest neutralino is in most models the hghtest 
supersymmetric particle (LSP) and usually just to be called as "the neutrahno" . Note 
here that, although the hghtest neutralino is the most studied LSP and some authors 
even use the LSP to indicate it, there are also some other candidates for LSP, e.g., the 
gravitino. 

Furthermore, a i?-parity should be also presented here [1]: 

i? = , (2.1) 

where B and L are the baryon and lepton number, S is the spin. For ordinary particles 
R = +1 (even) and for supersymmetric particles R = —1 (odd). If the i?-parity is 
conserved, one SUSY particle (with R = —1) could only decay to a lighter SUSY particle 
(with R = —1) and any number of ordinary particles (with R = +1). Certainly, such 
decay can not happen with the "lightest" SUSY particle since no SUSY particle can be 
hghter than the hghtest one. Hence, in models with strict i?-parity conservation, the LSP 
must be absolutely stable and is then the best candidate for Dark Matter [1]. 

In contrast, if it!-parity is broken, there is no special selection rule to prevent the 
decays of the supersymmetric particles in the spectrum with masses of order of a few 
GeV or larger [1]. In particular, there were no natural candidate for Cold Dark Matter 
[1]. Theories with broken i?-parity also possess baryon- and lepton- number violating 
interactions with strengths controlled by the scale of i?-parity violation [1]. 

2.1.2 Weakly Interacting Massive Particles (WIMPs) 

Weakly Interacting Massive Particles (WIMPs) x arise e.g., in supersymmetric exten- 
sions of the Standard Model of electroweak interactions and are the leading non-baryonic 
candidates for Cold Dark Matter [1] . They are stable particles and interact with ordinary 
matter only via weak interactions. Typically their masses have been presumed to be be- 
tween 10 GeV and a few TeV [1], [37]. They could include neutrahnos, sneutrinos, heavy 
fourth-generation Dirac and Majorana neutrinos, and non-minimal neutralinos (neutrali- 
nos in non- minimal supersymmetric models) [1]. 

Relic elementary particles are left over from the Big Bang. Stable or long-hved par- 
ticles with very weak interactions can remain in sufficient numbers to account for a 
significant fraction of critical density [15]. Very weak interactions are necessary for their 
annihilations to cease before their numbers are too small [15]. 

WIMPs exist in thermal equilibrium and in abundance in the early Universe, when the 
temperature of the Universe T exceeds their masses (T > m^) [1]. The equilibrium 
abundance is maintained by pair annihilation of WIMPs with their antiparticles x iiito 
lighter particles / (quarks and leptons, or even gauge- and Higgs-bosons if their masses 
are heavy enough) and vice versa {xX ^ ^0 [!]• The rate of this reaction is proportional 
to the product of the WIMP number density and the WIMP pair annihilation cross 
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(2.3) 



section into SM particles a a times the relative velocity between the two WIMPs in their 
center-of-mass system v [1] , [37] : 

Tx = n^{(^Av) , (2.2) 

where (■ ■ ■) indicates the thermal averaging. 

As the Universe cools to a temperature less than the masses of WIMPs (T < m^), 
the equilibrium abundance (number density) of WIMPs drops exponentially until the 
rate for the annihilation reaction {xx ~^ ^0 becomes smaller than the Hubble expansion 
rate of the Universe (F^ < H). At this point the interactions which maintain thermal 
equilibrium "freeze out" and the WIMPs cease to annihilate and drop out of thermal 
equilibrium [37]. Hence, a relic cosmological abundance "freezes in" [1], i.e., the density 
of the CO- moving WIMPs remains essentially constant. 

The time evolution of the number density of WIMPs ny^{t) can be described by the 
Boltzmann equation [1]: 

~dt 

Here is the number density of WIMPs in thermal equilibrium. The second term on the 
left-hand side accounts for the expansion of the Universe, the first term in the brackets on 
the right-hand side accounts for the depletion of WIMPs due to their pair-annihilation, 
and the second term arises from creation of WIMPs from the inverse reaction [1]. In 
the absence of number-changing interactions, the right-hand side would be zero and we 
would find n oc l/a^(t) [1], where a,{f) is the scale factor of the Universe in Eq.(1.12). 

Note that Eq.(2.3) describes both Dirac particles ^ as well as Majorana particles 
[1], which are particles that themselves are also their antiparticles, such as neutralinos 
(x = x) [!]• For the case of Majorana particles (they are so-called "self-annihilating"), 
the annihilation rate in Eq.(2.2) should be modified to [1] 

Tx = (^) M . (2.4) 

However, in each annihilation, two particles are removed and the factor of 2 can be 
canceled. For Dirac particles with no particle-antiparticle asymmetry, i.e., = n^, 
Eq.(2.2) is true. But the total number of particles plus antiparticles must be 2n^ [1]. In 
the case of Dirac particles with a particle-antiparticle asymmetry, the relic abundance is 
generally that given by the asymmetry [1] . For example, the relic proton density is fixed 
by the proton- antiproton asymmetry, i.e., the baryon number of the Universe [1]. 

The early Universe is radiation dominated and the Hubble expansion rate falls with 
temperature as [1] 

H{T) = 1.66 f . (2.5) 
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on the right-hand side should be modified to n-^^n^ for Dirac particles, but often n-^^ = n-^ has 
been assumed. 
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Here Mpi is the Planck mass [11] 



/ he 

Mpi = J— = 1.2209 X 10^^ GeV/c^ = 2.1764 x 10"^ g , (2.6) 
V G^N 

and the quantity is the effective number of relativistic degrees of freedom. It is approx- 
imately equal to the number of bosonic relativistic degrees of freedom plus | times the 
number of fermionic relativistic degrees of freedom [1] . While, for very high temperature 
(T>m^) [1], 

n'^ocT^. (2.7) 

Hence, the expansion rate H{T) in Eq.(2.5) decreases less rapidly than the number den- 
sity of WIMPs. This means that, at early times, the expansion term in Eq.(2.3), 3Hn^, 
is negligible compared with the right-hand side, and the number density tracks its equi- 
librium abundance [1]. 

However, at later times or at low temperatures (T < m^), the right-hand side in 
Eq.(2.3) becomes negligible compared with the expansion term, and the co-moving abun- 
dance of WIMPs remains unchanged [1]. It can be found that [1] 

^7-9[^) e-^l\ (2.8) 

where g is the number of internal degrees of freedom of the WIMPs and thus their density 
is Boltzmann suppressed [1]. If the expansion of the Universe were so slow that thermal 
equilibrium was always maintained, the number of WIMPs today would be exponentially 
suppressed (essentially, there would be no WIMPs) [1]. The temperature Tp at which the 
WIMPs freeze out is given by r^(Tp) = H{Tf) [1]. Using typical weak-scale numbers, 
the freeze-out temperature turns out to be [1] 

TV ~ ^ . (2.9) 
20 ^ ' 

There is a small logarithmic dependence on the mass and annihilation cross section here 
[1]. As stated above, after freeze out, the abundance of WIMPs per co-moving volume 
remains constant. 

Finally, the present relic density of WIMPs is then approximately given by (ignoring 
logarithmic corrections) [1], [37] 

^,2 0.1c pb 3xlO-2^cmVs 

n^h^ ~ const. X — 2 — - ~ ^ = ^ , (2.10 

where Tq is the current CMB temperature given in Eq.(1.22), and c is the speed of light. 
It is inversely proportional to the annihilation cross section of WIMPs. Hence, as the 
annihilation cross section is increased, the WIMPs stay in equilibrium longer, and we 
are left with a smaller relic abundance [1]. The annihilation cross section is generally 
expected to decrease as the WIMP mass is increased, so the relic abundance should be 
also increased [1]. Therefore, heavier WIMPs should be more likely to contribute too 
much to the mass of the Universe, and then be cosmologically inconsistent [1]. 
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2.1.3 Neutralinos 



As introduced in Subsec. 2.1.1, neutralinos are linear combinations of photino, Z-ino 
and neutral higgsinos (the supersymmetric partners of the photon, and neutral Higgs 
bosons, see Table 2.1): 

X° = aa + hZ + + diH'' , i = 1, 2, 3, 4. (2.11) 

In most SUSY models, the lightest neutralino is the lightest supersymmetric particle [37] 
and therefore the best motivated and the most widely studied candidate for WIMP Dark 
Matter, but not the unique candidate for LSP (e.g., sneutrinos) [37]. 

There are some theoretical reasons to believe that the lightest neutralino should be 
the LSP. First, suppose a charged uncolored SUSY particle, such as a chargino or a 
slepton, were the LSP. The relic number density of such particles can be given as roughly 
10"^nbM/GeV [1], where nb is the baryon number density and M is the mass of such 
particles. Then they would show up in searches for anomalously heavy protons [1]. 
Null results from such searches rule out such charged particles over a broad mass range 
[1]. Moreover, grand unified models predict relations between the masses of the SUSY 
particles . In most models the gluino is more massive than the neutralino, and the squarks 
are also heavier than the sleptons [1]. In addition, some detailed calculations show that 
the lightest neutralino has the desired thermal relic density, I^dm in Ect-(l-47), in at least 
four distinct regions of parameter space [37] . 

2.1.4 Sneutrinos 

Sneutrinos are the spin-0 supersymmetric partner of the neutrinos (see Table 2.1). 
There are some reasons to rule out sneutrinos to be good candidate for Dark Matter. 
First, in most models, there is a slepton with mass similar to, but slightly smaller than, 
the sneutrino mass [1]. Meanwhile, their masses would have to exceed several hundred 
GeV for them to make good Dark Matter candidates. This is uncomfortably heavy for the 
lightest sparticle [37]. On the other hand, the annihilation cross sections of sneutrinos 
arc expected to be quite large [37]. Hence, the negative outcome of various WIMP 
searches rules out ordinary sneutrinos as primary component of the Dark Matter halo of 
our Galaxy [37]. However, in models with gauge-mediated SUSY breaking the lightest 
messenger sneutrino could be a good candidate [37]. 

2.1.5 Heavy fourth-generation Dirac and Majorana neutrinos 

They are the first proposed WIMP candidates for CDM [1] . They are heavy, but stable 
particles and assumed to have (weak) interactions with ordinary matter though Standard 
Model coupling to the Z^ boson [1]. Such neutrinos could annihilate into light fermions 
via s-channel exchange of a Z^ boson. For m ^ Mz the cross section is proportional to 
the square of their mass [1]. Because their interactions are fixed by gauge symmetry, the 
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only adjustable scale is then their masses [1]. The cosmological abundances of the heavy 
Dirac and Majorana neutrinos have been given as [1] 



for neutrino masses in the range 0{1 GeV) < rrii, <^ niz = 91.19 GeV. 

However, there is no obvious reason why such massive neutrinos should not be allowed 
to decay [37]. Moreover, an SU(2) doublet neutrino will have a too small rehc density if 
its mass exceeds as required by LEP data [37]. 

On the other hand, for such neutrinos with masses greater than the electroweak gauge- 
boson masses, annihilations into gauge- and/or Higgs-boson pairs could occur. However, 
the cross section would not decrease as the neutrino mass increases [1], so the relic 
abundance of neutrinos with masses of the order of 100 GeV remains too small to account 
for the Dark Matter in the Galactic halo [1] . 

Dirac neutrinos interact with nuclei through a coherent vector interaction [1] (some 
details about the vector interaction with nuclei will be given in Subsec. 3.4.1). Thus 
the Dirac- neutrino- nucleus cross section is expected to be quite substantial [1], and this 
would lead to a significant event rate in a direct detection experiment. Null results from 
such experiments have ruled out Dirac neutrinos with masses in the range 12 GeV < m,y^D 
< 1.4 TeV as the primary component of the Dark Matter halo [1]. 

Meanwhile, Majorana neutrinos interact with nuclei only via an axial-vector inter- 
action [1] (some details about the axial-vector interaction with nuclei will be given in 
Subsec. 3.4.2), and are therefore difficult to detect directly. However, such neutrinos 
would be captured in the Sun by scattering from hydrogen therein and their pair anni- 
hilations in the Sun would produce energetic neutrinos from the Sun [1]. Null results 
from searches for energetic neutrinos at e.g., Kamiokande have also ruled out Majorana 
neutrinos with mass less than a few hundred GeV [1] . 

2.1.6 Axions 

Axion a is also one of the leading candidates for CDM. Axions have been introduced by 
Peccei and Quinn to solve the strong CP (charge-conjugation and parity) violation prob- 
lem of QCD. They are pseudo Nambu-Goldstone bosons associated with the spontaneous 
breaking of a new global Peccei-Quinn (PQ) U(l) symmetry at scale /„ [37]. 
The present relic density of the axions can be given as [37] 




(2.12) 



and 




(2.13) 




(2.14) 
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where is a numerical factor lying roughly between 0.5 and a few, 9a is a "misalignment 
angle" which parameterizes the axion field. Suppose 9a ~ C'(l), axions will have the 
required cosmological energy density in Eq.(1.47) to be Dark Matter, if /„ ~ 0(10^^ GeV). 
It is pretty comfortably above the laboratory and astrophysical constraints and this would 

correspond to an axion mass iria ~ 10~^ eV [37]. 

Axions could be detected by looking for their conversion to microwave photons, a — > 7, 
in a strong magnetic field [37]. Such a conversion could proceed through the loop-induced 
077 coupling, whose strength Qaj^y is thus an important parameter of axion models [37]. 
Moreover, the conversion rate can be enhanced in a high quality cavity on resonance and, 
due to the equation rriaC^ = hcUres, varying this resonance frequency can give a range of 
rua, or, equivalently, /„ [37]. 

2.1.7 Other possible SUSY candidates 

Besides the neutralinos and the sneutrinos, some other supersymmetric particles are 
also (theoretically) possible to be candidates for Dark Matter. ^ 

Axino is the spin-| superpartner of the axion. It may be the LSP or the next-lightest 
supersymmetric particle (NLSP) and may decay to the LSP [1]. When the axino is the 
lightest supersymmetric particle and has a mass of a few keV, it can be a good candidate 
for Warm Dark Matter (WDM) [1]. 

While, gravitino, the spin-| superpartner of graviton, is also a possible candidate for 
Dark Matter. The gravitinos will decouple at temperatures of order of the Planck scale 
{Epi ~ 10^^ GeV [1]). Thus the physics of the gravitinos must be considered at energies 
and temperatures right up to this scale [1]. In addition, if gravitinos behave as standard 
stable thermal relics with an abundance determined by consideration of their decoupling, 
the mass of gravitinos should be less than a few keV [1]. However, in some models with 
gravitinos as LSP, the NLSP should decay to a gravitino plus ordinary particles [1] . Since 
the coupling to gravitinos is so weak, this NLSP will be very long-lived and the products 
of its decay will contain 7-ray with high energies [1] . 

2.2 Hot Dark Matter (HDM) 

"Hot" has been used here to indicate that such Dark Matter particles moved relativis- 
tically in the early Universe. Due to their fa,st velocities, they would cover great distances 
and then form some very large scale structures. This means that the Hot Dark Matter 
forms the structure of our Universe from the top down, with superclusters fragmenting 
into clusters and galaxies [15]. It is in contrast to the observational evidence which indi- 

^Bcsidcs the different supersymmetric extensions of the Standard Model, there are also some theories 
based on "flat universal extra dimensions (UED)" . The most studied candidate for CDM in these extra- 
dimension models is the first Kaluza-Klein (KK) mode of the hypercharge gauge boson (the lightest KK 
particle, LKP) 7(1). 
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cates that the structure of our Universe has been formed from the bottom up by merging 
dust to galactic scale structures [15]. 

However, there are still some suggestions in which part of Dark Matter is hot and the 
rest is cold. In these models the bulk of the Dark Matter (especially in galactic halos) is 
still cold. 

2.2.1 Massive neutrinos 

The leading candidates for Hot Dark Matter are the massive neutrinos. As shown in 
Subsec. 1.2.9, WMAP results combined with other astronomical measurements lead to a 
contribution for light (but massive) neutrino species [11]: 

0.014. (1.51) 

They could include the electron-, muon-, and tauon-neutrinos in the Standard Model 
with non-zero masses ^ as well as the forth-generation Dirac and Majorana neutrinos 
(described in Subsec. 2.1.5) with extremely light masses. 

2.3 Dark baryons 

As mentioned above, some CDM particles, e.g., neutralinos and axions, could form 
galactic scale structures such as galaxies and clusters of galaxies, while, some RDM 
particles, e.g., massive neutrinos, could form larger structures of the Universe. This 
means that on different scales Dark Matter might consist of different materials [1]. ^ 

Moreover, in this chapter I presented some theoretically predicted (SUSY) particles 
as candidates for Dark Matter. However, until now there is no direct accelerator evidence 
for the existence of supersymmetry [1]. Actually, it is not absolutely certain that Dark 
Matter is neither baryons nor neutrinos [1]. There are also some conservative cosmological 
models which describe the Universe only in terms of baryons and perhaps neutrinos [1]. 

On the other hand, as shown in Subsec. 1.2.9, the baryonic matter density in the 
Universe is 

Ob - 0.042, (1.46') 
but only around 25% of the baryonic matter are luminous: 

OlumC^O.Ol. (1.48) 

Although, as mentioned in Subsec. 1.2.6, most of the baryons in the clusters of galaxies 
reside not in the galaxies themselves, but in form of hot intercluster, x-ray emitting gas 



^At present we know from v oscillations that at least two of these three SM neutrinos have small, but 
non- vanishing masses. 

^However, some recent researches indicate that there should be only one species of Dark Matter in 
the Universe. 
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[15], such hot gas in the clusters of galaxies only accounts for around 10% of the baryons 
in the Universe [15]. Hence, there should (must) be some baryonic Dark Matter. ^ 

Two most promising possibilities for such dark baryons are diffuse hot gas and dark 
stars, which include white dwarfs, neutron stars, black holes, or objects with masses 
around or below the hydrogen- burning limit [15]. 

2.3.1 Massive astrophysical compact halo objects (MACHOs) 

Massive astrophysical compact halo objects include, for example, brown dwarfs which 
are balls of hydrogen and hehum with masses < O.OSMq and therefore never begin nuclear 
fusion of hydrogen [1] (but they do burn deuterium), jupiters which are similar to brown 
dwarfs but have masses ~ O.OOIM© [1] and do not burn anything, and white dwarfs 
[1]. Actually, objects with masses around or below the hydrogen-burning limit could be 
baryonic Dark Matter [15]. 

Meanwhile, neutron stars and stellar black-hole remnants are also candidates for bary- 
onic Dark Matter [1]. Black holes with masses ~ IOOMq could be remnants of an early 
generation of stars which were massive enough so that not many heavy elements were 
dispersed when they underwent their supernova explosions [1]. Primordial black holes 
which formed before the era of Big Bang could be counted for non-baryonic Dark Matter 
rather than baryonic one [37]. However, such an early creation of a large number of black 
holes is possible only in certain somewhat contrived cosmological models [37]. 

MACHOs might represent a large part of the galactic Dark Matter and could be 
detected through the microlensing effect [37]. The MACHO, EROS, OGLE collaborations 
have performed programs of observation of such objects by monitoring the luminosity of 
millions of stars in the Large and Small Magellanic Clouds [37]. They concluded that 
MACHOs contribute < 40% (MACHO) or even < 20% (EROS) to the mass of the galactic 
halo [37]. 



'"'However, some recent results of the measurement of the opacity of the Lyman-a forest toward high- 
redshift quasars indicate that there are probably enough baryons at > 3, but it is not clear where they 
are now. 
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Chapter 3 

Direct Detection of WIMPs 



By definition, Dark Matter could neither emit nor absorb electromagnetic radiation. 
However, as described in Subsec. 2.1.2, WIMPs would have annihilated to some ordinary 
matter, e.g., quarks and leptons, in the early Universe. Otherwise, they would have 
unacceptable large abundance today. According to the crossing symmetry, the amplitude 
for WIMP annihilation to, for example, quarks is related to the amphtude for elastic 
scattering of WIMPs from quarks [1]. Therefore, WIMPs should have some small, but 
non-zero couplings to ordinary matter. 

Due to this couphng to nucleus (through the coupling to quarks), WIMPs could 
scatter elastically from target nuclei of the detector material and produce nuclear recoils 
which deposit energy in the detector. Hence, one of the most promising methods of 
detecting Galactic Dark Matter is the direct detection of WIMPs [38]-[44]. Note that, 
although the lightest neutralino is the leading candidate for Dark Matter, such WIMP 
direct searches are not specialized to detect the neutralino but any particle with similar 
generic properties, e.g., a mass between a few GeV and a few TeV and weakly interacting 
with ordinary matter [45]. ^ 



3.1 Elastic WIMP- nucleus scattering 

Using the standard assumption of the WIMP density near the Earth 

po ^ 0.3 GeV/cVcm^ , (1-55) 
and assuming a WIMP mass 

m^^lOO GeV/c^ , (3.1) 

^The lightest Kaluza-Klein particle arising in the extra-dimension models (mentioned as footnote in 
Subscc. 2.1.7) can also scatter elastically from the detector miclci through KK-quark qi^i^ and Higgs 
exchange [46]. Thus such particles could also be detected from direct detection experiments. A brief 
description about the interaction between (/(i) and Higgs and the analysis using recent experimental 
results can be found in Ref. [46]. 
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the number density of WIMP can be found as 



n 



a; 3 X 10"^ cm 



-3 



(3.2) 



Meanwhile, by the assumption that the halo WIMPs are gravitationally bound to the 
Galaxy and its halo, the average velocity of WIMP wind is then approximately equal to 
the stellar velocity in the Solar neighborhood: ^ 



Therefore, the WIMP flux is ~ 10^ WIMPs per square centimeter of the Earth's surface 
per second. 

However, the very low cross section of WIMPs on ordinary material makes the elastic 
WIMP-nucleus scattering very rare. In typical SUSY models with neutralino WIMPs, 
WIMP-nucleus cross section is about 10"^ ~ 10"^ pb (lO"''^ ~ 10"''° cm^) ^ and the 
expected event rate is then at most 1 event/kg /day [37], in some models it is even less 
than 1 event/ton/yr [47]. 

With expected WIMP mass in the range 10 GeV/c^ to 10 TeV/c^ [37], typical nu- 
clear recoil energies are of order of 1 to 100 keV. However, as we can see in Fig. 4.1 in 
Subsec. 4.2.1, the event rate drops approximately exponentially and most events should 
be with energies less than 40 keV (a simple theoretical estimate will be given in Sub- 
sec. 3.5.1). 

On the other hand, in the energy range from a few to a couple hundred keV, typical 
background noise due to cosmic rays and ambient radioactivity is much larger. Thus a 
underground laboratory and extensive shielding around the detector to protect against 
cosmic-ray induced backgrounds, and selection of extremely radiopure materials are nec- 
essary and important [37] (more details about background and its discrimination will be 
given in Sec. 3.6). 

The event rate of elastic WIMP-nucleus scattering depends on various parameters 
coming from astrophysics, particle physics and nuclear physics: the WIMP density near 
the Solar system po, the WIMP-nucleus cross section, the WIMP mass m^, and the 
velocity distribution of the incident WIMPs f{v) in the Galactic halo near the Earth. 
However, by some standard assumptions about the halo model, e.g., the WIMP density 
profiles (have been presented in Sec. 1.3) and velocity distributions (will be discussed in 
Subsecs. 3.1.3 and 3.2.1), the expected event rate mainly depends on two unknowns: the 
mass of the incident WIMPs and the WIMP-nucleus cross section. Hence the experimen- 
tal observable is usually expressed as a contour in the WIMP mass-cross section plane 
(e.g., Figs. 3.4 and 3.5), although it is basically the scattering rate and is a function of 
energy [37]. 

^Notc that (v) = (|v|) is 0{vq) even though (v) = (see footnote on p. 15). 
31 barn = IQ-^^ crn^, 1 pb = IQ-^e cm^ 



(v) fti 250 km/s . 



(3.3) 
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3.1.1 Rate for elastic WIMP-nucleus scattering 

The direct detection experiment measures the number of events per unit time per 
unit mass of detector material as a function of the energy deposited in the detector Q. 
Quahtatively, the event rate of direct detection, i?, can be simply expressed as [1] 

^ n(v) a ,„ 
R ^ ^ ' , (3.4 

where {v) is the average velocity of the incident WIMPs relative to the Earth frame (i.e., 
to the target), a is the WIMP-nucleus cross section, and is the mass of the target 
nucleus. Here we multiply the factor l/m^ to get the number of target nuclei per unit 
mass of the detector material. 

More accurately, one should take into account the fact that the WIMPs move in the 
halo with velocities determined by their velocity distribution function f{v), and that the 
differential cross section depends on f{v) through an elastic nuclear form factor F{q) [1]: 

Here ctq is the total cross section ignoring the form factor suppression, 

rrij. — (3.6) 

m-)^ + mN 

is the reduced mass, q is the transferred 3-momentum: 

q = pm^Q . (3.7) 

Therefore, in general, the differential scattering event rate (per unit detector mass) should 
be written as [1] 

dR = — — — / vfi (v) da dv 

'fiiv) 



W)/ 



dvdQ, (3.8) 



where fi (v) is the one-dimensional velocity distribution function of WIMPs impinging on 
the detector, v is the absolute value of the WIMP velocity in the Earth rest frame, and we 
have to integrate over all possible incoming velocities. By means of classical mechanics, 
the transferred momentum q can be expressed as 



q^2 



mN 



sin(gf^)^2m.J ^-7^^" , (3.9) 



rriy^ + mN 

where ^cm is the scattering angle in the center-of-momentum frame. Since 
< 1 -COS^CM < 2, 
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for a given deposited energy Q, we have 

2mN ?7T.N ' 

i.e., the minimal incoming velocity of incident WIMPs that can deposit the energy Q in 
the detector can be expressed as 



M-J^,^/Q-»y/Q, (3.10) 



where I have defined 



Then the differential event rate for elastic WIMP-nucleus scattering, Eq.(3.8), can be 
rewritten as 



dR 
dQ 



AF\Q) / 

Jv„ 



dv, (3.12) 



where the constant coefficient A is defined as 

A^:^. (3.13) 

Note that, first, a defined in Eq.(3.11) depends only on the WIMP mass (and the 
mass of the target material, m^, which we can choose). The two as yet unknown pa- 
rameters, i.e., the WIMP density po and the total WIMP-nucleus cross section (Jq, have 
been collected in the coefficient A defined in Eq.(3.13). Second, I assumed here that the 
detector essentially only consists of nuclei of a single isotope. If the detector contains 
several different nuclei (e.g., Nal as in the DAMA detector [48]), the right-hand side of 
Eq.(3.12) has to be replaced by a sum of terms, each term describing the contribution of 
one isotope. For simplicity, in the remainder of this work I will focus on mono-isotopic 
detectors. 

Finally, the total event rate per unit time per unit mass of detector material can be 
expressed as 

where Qthre is the threshold energy of the detector. 



3.1.2 Nuclear form factor (for spin- independent coupling) 

Here I present two most commonly used parameterizations of the squared nuclear 
form factor F'^{Q) in Eq.(3.12) for spin-independent coupling, which usually dominates 
the event rate (more details about WIMP-nucleus couplings will be given in Sec. 3.4). 
Moreover, the form factors for spin-dependent coupling are still only poorly understood. 
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The simplest form factor is the exponential one, first introduced by Ahlen et al. [49] 
and Preese et al. [42]: 



,-Q/Qo 



(3.15) 



where Q is the recoil energy transferred from the incident WIMP to the target nucleus, 
1.5 



Qo 



is the nuclear coherence energy and 



Rq — 



0.3 + 0.91 



VGeV; 



fm 



(3.16a) 



(3.16b) 



is the radius of the nucleus. The exponential form factor implies that the radial density 
profile of the nucleus has a Gaussian form. This Gaussian density profile is simple, but 
not very realistic. Engel has therefore suggested a more accurate form factor [50], inspired 
by the Woods-Saxon nuclear density profile, 

,n 2 



3ji(gi?i 

qRi 



-{is? 



(3.17) 



Here ji{x) is a spherical Bessel function, q is the transferred 3-momentum given in 
Eq.(3.7), and 



R\ - 5s2 



with 

Ra ^ 1.2 A^'^ im, s~lfm, 
where A is the atomic mass number of the nucleus. 



(3.18a) 



(3.18b) 



3.1.3 Simple isothermal Msixwellian halo 

For the simplest halo model presented in Subsec. 1.3.2, the canonical isothermal 
spherical halo, with the assumption that the WIMPs trapped in the galactic field have 
attained thermal equilibrium with a Maxwellian velocity distribution [51], the velocity 
distribution function is given by [1] 



(3.19) 



where is the orbital velocity of the Sun in the Galactic frame: 

vo~ 220 km/s, (1.53') 

which characterizes the velocity of all virialized objects in the Solar vicinity. Then, since 

Sv = v^dv dVt = Anv'^dv , 
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the normalized one-dimensional velocity distribution function has been obtained as [1] 



/.,o..W-;|(9e-«V< (3.20) 

According to Eq.(3.12), the scattering spectrum of the simplest theoretical velocity dis- 
tribution given in Eq.(3.20) can be obtained as (a detailed calculation will be given in 
App. C.1.1) 

Meanwhile, the mean velocity and velocity dispersion of the halo WIMPs can be obtained 
as (detailed calculations will be given in App. B.1.1), respectively, 

■ 2 ■ 



(^^)Gau = ^ 'y/l,Gau(t') dv = f — = j , (3.22) 



and 

(^^')Gau = f v'fi,GUv) dv = Q v', . (3.23) 

For light WIMPs, the effect due to the form factor introduced in Eq.(3.5) can be neglected 
and we can use F'^{Q) ^ 1 [1]. Then the total event rate in Eq.(3.14) can be found directly 
as 

i?Gau(Qth.e) = ^ (^] e-^^^^-^/^o . (3.24) 

For the case of Qthre — 0, this result can be reduced to 

i?Gau(gth.e = 0) = ^2^^ , (3.25) 

which is exactly the naive estimate in Eq.(3.4). On the other hand, with the exponential 
form factor F^^{Q) given in Eq.(3.15), one can find that (a detailed calculation will be 
given in App. C.2.1) 

i?Gau,e.(gth.e) = ^^^^ ^-'^^Q.^.J^l^^) , (3.26) 

and then 

-RGau,ex(Qthre = 0) = ^0 o( )Gau _ ^2 ^ (3.27) 

where I have defined 

/ 2 \ -1/2 

It can be seen that, for the case that the exponential form factor F^-^{Q) can be neglected, 
or, equivalently, Qo — oo, i.e., /5 — >• 1, -RGau,ex(Qthre) in Eqs.(3.26) and (3.27) will reduce 
to Eqs.(3.24) and (3.25). 
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WIMP 



Wind 




Figure 3.1: The Earth orbit in Galactic coordinates. The Sun moves to the left with 
about 220 km/s, inducing a WIMP wind (figure from [45]). 

3.2 Annual modulation of the event rate 

The one-dimcnsional velocity distribution function /i,Gau('^) given in Eq.(3.20) has 
been considered in the Galactic rest frame. More realistically, the orbital motion of the 
Solar system in the Galaxy as well as the motion of the Earth around the Sun must be 
considered [41]- [43]. As shown in Fig. 3.1, since the speed of the Earth adds to or subtracts 
from the speed of the Sun, the event rate for a given recoil energy (or energy range) in 
Eq.(3.12) should be a cosinusoidal function with a one-year period (see Figs. 3.3 and 5.1) 
and a peak around June 2nd [42] . The expected amplitude of this annual modulation is 
around 5%. ^ 

Originally, such an annual modulation was only expected for the WIMP signal, not for 
the background. Thus this effect might serve a method to distinguish the WIMP signal 
from the background. And actually, the DAMA collaboration [48] has claimed that they 
have observed this annual modulation of the event rate [52]- [54] (more details about the 
DAMA result will be given in Subsec. 3.7.3). However, the much larger background might 
also be subject to modulation [37]. For example, the dependence of the cosmic muon fiux 
on the atmospherical temperature, or the dependence of the background neutron flux on 
water in rock and concrete. Hence, the signal identification should also be performed. 

"^Thc ratio of a theoretically expected amplitude of this annual modulation to the time-averaged 
scattering spectrum as function of the recoil energy will be given in Fig. 5.2 in Subsec. 5.2.1. It can be 
seen that, for recoil energy between and 50 keV, the modulated amplitude is around —4% 5%. Some 
detailed discussions about the annual modulation of event rate will be given in Sees. 5.1 and 5.2.1. 
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3.2.1 Shifted Maxwellian halo 



When we take into account the orbital motion of the Solar system around the Galaxy, 
as well as that of the Earth around the Sun, the velocity distribution function in Eq.(3.20) 
should be modified to [1] 



fl,sh{'^,Ve) = 



with 



TT WeVo 



1.05 + 0.07 cos 



'2'K{t - tp 



1 yr 



(3.29) 



(3.30) 



where tp ~ June 2nd is the date on which the velocity of the Earth relative to the WIMP 
halo is maximal [42]. Eq.(3.30) includes the effect of the rotation of the Earth around 
the Sun (second term), but does not allow for the possibility that the halo itself might 

rotate (some discussions about such bulk rotation has been given in Subsec. 1.3.6). 

Substituting the shifted Maxwellian velocity distribution function in Eq.(3.29) into 
Eq.(3.12), the theoretically expected scattering spectrum can be obtained as (a detailed 
calculation will be given in App. C.1.2) 



sh 



erf( "^+^- ) -erf( "^-^- ) 



(3.31) 



Here erf(x) is the error function, defined as 
2 



erf(x) = [ e *^dt. 

VTT JO 



Meanwhile, the mean velocity and velocity dispersion in Eqs.(3.22) and (3.23) should be 
modified to (detailed calculations will be given in App. B.1.2), respectively. 



{v) 



sh 



TT 



and 



2 I 2 



(3.32) 



(3.33) 



As what I did in Subsec. 3.1.3, considering the hght-WIMP case and using F'^{Q) f» 1, 
the total event rate for the shifted MaxweUian distribution in Eq.(3.29) can be found as 
(a detailed calculation will be given in App. C.1.2) 



thre I 



myrriTsi \2ve 



cd{S+) - erf(5_) 



TT 



(3.34) 
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where I have defined 

r, _ aVQthre ± V, 
^± = 



o± = . 

For the case of Qthre — 0, -Rsh(Qthre) in Eq.(3.34) can be reduced directly 

/tt, 



(3.35) 



to 



-Rsh(<5thre — 0) — 

where I have used 

erf(— x) = — erf(a;) . 



(3.36) 



Moreover, for the case with the exponential form factor F^^{Q) gi 
(a detailed calculation will be given in App. C.2.2) 



;iven in Eq.(3.15), I have 



-'•'shjexVVthre j — 



V, 



'0 



f3' 



myjUfi \2veJ \l — P"^ J 



erf (5+) - erf (5_) 



- /?e-(^-^')^'/''o [erf (T+) - erf (T.; 



where I have defined 



, (3.37) 



(3.38) 



For the case of Qthre = 0, i?sh,ex(Qthre) in Eq.(3.37) can be reduce to 



-fi'sh,ex(Q 



thrc 



0) 



erff^)-/5e-M^K/.^erf (^] 



m^m^ \VeJ \i — p ' * ' \ ; 

It is not difficult to check that -Rsh,ex(Qthre) in Eqs.(3.37) and (3.39) can be reduced to 
Eqs.(3.34) and (3.36) when one neglects the form factor F^J^Q), i.e., let (3^1. On the 
other hand, as v^. <^ fg, one can also prove that the results in Eqs.(3.32) to (3.34), and 
(3.36) can be reduced to Eqs.(3.22) to (3.25). 



3.3 Diurnal modulation of the event rate 

Similar to the annual modulation caused by the orbital motion of the Earth around 
the Sun, due to the rotation of the Earth, the event rate for a given energy (or energy 
range) should have a diurnal modulation [45], [22]. 

There are two different effects caused by this diurnal modulation. The first one is 
the shielding of the detector of the incident WIMP flux by the Earth [22] (illustrated 
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(b) 



Figure 3.2: Two effects caused by the diurnal modulation: (a) shielding of the detector 
of the incident WIMP flux by the Earth (figure from [22]), and (b) directionality of the 
WIMP wind (figure from [45]). 
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in Figs. 3.2(a)). Some authors claimed that, for WIMP masses close to 50 GeV/c^ and 
under certain assumptions, this diurnal modulation due to the shielding of the WIMP 
flux could be larger than the annual modulation [43], [22]. However, this requires a large 
WIMP-nucleus cross section and recent experimental results have (almost) excluded this 
possibility. Moreover, practically, it is still impossible for the detectors nowaday and 
should also be very difficult for the next-generation ones to get more than a few singles 
per day to prove this effect (more details about the status of the operated experiments 
and their results will be given in Sees. 3.7 to 3.9). 

On the other hand, the second effect due to the rotation of the Earth is the direc- 
tionality of the WIMP wind: a daily forward/backward asymmetry of the nuclear recoil 
direction (illustrated in Figs. 3.2(b)). A gaseous detector (e.g., DRIFT [55]) or anisotropic 
response scintillators should have the ability to measure this recoil direction [45], [37]. 



3.4 Target material dependence 

The WIMP-nucleus cross section ao in Eqs.(3.5) and (3.13) depends on the nature 
of the WIMP couphngs to nucleons. For non-relativistic WIMPs, one in general has to 
distinguish spin-independent (SI) and spin-dependent (SD) couplings [37]. 

3.4.1 Spin- independent (SI) cross section 

The total cross section for "scalar" coupling can be expressed as [1] 



4m: 



2 



(T0,scalar = ^k/p + (A-^)/nl • (3.40) 

Here mr is the reduced mass of the WIMP and the target nucleus in Eq.(3.6), Z is the 
atomic number, i.e., the number of protons, A ~ Z is then the number of neutrons, /p 
and /n are the effective couplings of WIMPs to protons and neutrons, respectively. 

Here we have to sum over the couplings to each nuclcon before squaring because the 
wavelength associated with the momentum transfer is comparable to or larger than the 
size of the nucleus [51], the so-called "coherence effect". In most cases, the couplings to 
protons and neutrons are approximately equal [1], 

/n ^ /p . (3.41) 

Then the cross section for scalar interaction in Eq.(3.40) can be reduced to 

(70,scalar OC A'^ . (3.42) 

This means that, due to the coherence effect with the entire nucleus, the cross section 
for scalar interaction scales approximately as the square of the atomic mass of the target 
nucleus. Hence, higher mass nuclei, e.g., Ge or Xe, are preferred for the search for the 
scalar interaction [37]. 
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On the other hand, WIMPs could also have a "vector" couphng to protons and neu- 
trons [1]: 

,2 



"^0, vector 



^^^^2Zb^ + {A-Z)bn\ , (3.43) 

where bp and b^ are the effective couplings to protons and neutrons. However, for Majo- 
rana WIMPs (x = x), e.g., the neutralino, there is no such vector interaction [37]. 

3.4.2 Spin-dependent (SD) cross section 

WIMPs could also couple to the spin of the target nucleus, an "axial-vector" (spin-spin) 
interaction. For this spin-spin coupling, only unpaired nucleons contribute significantly to 
the interaction, as the spins of the A nucleons in a nucleus are systematically anti-aligned 
[51]. And it is obvious that this spin-dependent interaction exists only if the incident 
WIMPs carry spin [37]. 

The total cross section for spin coupling can be expressed as [1] 

_ 32m2 r 

•^Ojaxial 



AV(J + 1) , (3.44) 

where J is the total angular momentum of the nucleus and A (oc 1/J) depends on the 
axial couplings of WIMPs to the quarks. 

Because of the dependence on the nuclear spin factor, the useful target nuclei for 
search for spin interaction are -'^F and ^^"''I [37]. 

3.4.3 Comparison of the SI and SD cross sections 

Generally speaking, a WIMP could have both scalar and spin-dependent interactions 
with the nucleus. Thus the WIMP-nucleus cross section (Jq in Eqs.(3.5) and (3.13) should 
be the sum of the scalar cross section, cro,scaiar) in Eq.(3.40) and the spin cross section, 
(7o,axiai, in Eq.(3.44). 

For the scalar interaction, an analytic nuclear form factor, e.g., the exponential and 
the Woods-Saxon form factors presented in Subsec. 3.1.2, can be used. For the spin 
interaction, the form factor will differ from nucleus to nucleus and no simple analytic 
form factor can provide a very good approximation [1]. Fortunately, for nuclei with 
A> 30 , the scalar interaction almost always dominates the spin interaction [1]. 



3.4.4 Target mass 

The scattering event rate depends also on the atomic mass of the target material 
directly. 

First, according to Eq.(3.10), the smaller a the lower the incoming velocity with which 
the incident WIMPs can deposit energy larger than the threshold energy. Meanwhile, 
according to the definitions of a and mr in Eqs.(3.11) and (3.6), it can be found that, for 
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WIMPs with a given mass and detector with a given threshold energy, a will be smallest 
if the mass of the target nucleus is equal to the WIMP mass m^. 

Second, for a given total mass of detector material, a larger target mass means also a 
smaller number density of the nucleus which can interact with the incident WIMPs. It 
will certainly reduce the total event number. 

3.5 Measurement of recoil energy 

3.5.1 A simple estimate 

As an example, we assume a WIMP mass 

m;^ fti 100 GeV/c^ (3.1) 

and use the standard theoretical WIMP rms velocity 

^ 270 km/s. (1.52) 

then the average kinetic energy of the incident WIMPs can be estimated as 

l,E^ ^ 40 keV . (3.45) 

On the other hand, by means of classical mechanics, the recoil energy of the target nucleus 
due to the elastic scattering can be expressed as 



Q 



Am^rriTss 2 



cos 6'Lab 



, (3.46) 



where ^^Lab is the recoil angle in the laboratory frame. This expression shows that the 
maximum recoil energy is obtained when = m^. This is also why this search should 
be more efficient for target material with a mass comparable to the WIMP mass. 

3.5.2 Induced signals 

When a WIMP scatters off a nucleus, the nucleus will at first obtain a few tens of kcV 
kinetic energy and then dissipate this energy in the detector via three main processes: the 
electrons can be stripped by the scattered nucleus and an ionized nucleus-electrons system 
will be produced, this electronic activity can emit light, and the movement of the recoiling 
nucleus in the lattice can also induce vibrational phonons. Moreover, the ionization and 
scintillation energy will convert into phonons that will eventually thermalize and produce 
a tiny elevation of the temperature in the detector. 

Hence, generally speaking, due to the elastic WIMP-nucleus scattering, the nuclear 
recoil can induce three different signals: ionization (charges), scintillation (light), and 
heat (phonons). 
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3.5.3 Quenching factor 

When a photon with energy between keV and MeV enters a detector, it will induce 
an electron recoil with a range of the order of the /xm and transfer most of its energy to 
the electron. However, the range of a nuclear recoil is only of the order of the nm and 
the nucleus will lose a substantial part of its energy directly into phonons associated with 
atom vibrations as the nucleus is stopped in the lattice [51]. 

Hence, the quenching factor (the nuclear recoil relative efficiency) for the ionization 
detectors has been defined as the ratio of the number of charge carriers produced by a 
nuclear recoil due to the WIMP interaction to that produced by an electron recoil with 
the same kinetic energy (energy calibrated with a 7-source, called "electron equivalent 
energy" or "eee"). Meanwhile, for scintillating detectors, the quenching factor is defined 
as the ratio between the light produced by a nuclear recoil and by an electron recoil. 
For conventional detectors, this factor is usually lower than 0.3 [22], [51]: ~ 0.3 for Ge 
or Si, ~ 0.25 for Na, ~ 0.09 for I, and ~ 0.2 for Xe. While, for cryogenic detectors 
measuring heat, the quenching factor has been measured to be around one for recoiling 
nuclei independently of the energy [22]. 

Note that, due to this quenching factor, the measured recoil energies are often quoted 
practically in keVee instead of true recoil energies in unit of keV. 

3.5.4 Heat 

Basically a cryogenic detector has been made of a crystal with a thermometer glued 
on it, and operated at very low temperature (around 20 mK). 

When the detectors have been cooled to the operating temperature in a dilution 
refrigerator, the heat capacity (oc T^) is so low that even a few keV of deposited energy 
raises the temperature of one of the detectors by a measurable amount, allowing the 
amount of energy deposited to be determined [1]. 

Moreover, a superconducting-normal phase transition due to the elevation of the tem- 
perature has been used by the CRESST collaboration [56]. A thin film of tungsten (W) 
can be grown on a silicon detector and held just below the critical temperature. Phonons 
created by a WIMP-nucleus scattering would heat the superconducting film, causing it 
to go normal, and the change in resistance could be measured [1] . A very low threshold 
energy (~ 500 eV) of such detector were reached by the CRESST-I experiment with a 
262 g sapphire detector [22], [57] (more details about the CRESST experiments and their 
results will be given in Subsec. 3.7.2). 

Similarly, it is also possible to use some small superconducting granules in a magnetic 
field as detector, when one of such detectors is heated by a nuclear recoil, it would go 
normal and thereby cause a measurable change in the magnetic flux [1]. 
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3.5.5 Ionization 

A small voltage is placed across the crystal of the detector, and when several atoms 
have been ionized, the freed electrons will drift to one side, the collected charges can be 
used as a measure of the energy deposited in ionization [1] . 

Germanium ^^Ge used initially in the neutrino-less double- /9 (Oz/2/9) decay experiments 
has been used as the first detector material for direct WIMP detection experiments by 
the Heidelberg-Moscow (HDMS) collaboration [58] (more details about the HDMS ex- 
periments and their results will be given in Subsec. 3.7.5). Thanks to the high intrinsic 
purity achieved by the semiconductor industries and the technique developed for the 
decay experiments, Ge ionization detectors have nowaday very low thresholds and very 
good resolutions (Qthre ^ 4: 10 keVge, equivalent to ~ 15 ~ 30 keV recoil energy, for 
HDMS [22]). Moreover, silicon (Si) has also been used by e.g., CDMS collaboration [59] 
as detector material (more details about the CDMS experiments and their results will be 
given in Subsec. 3.7.1). However, the size of such ionization detectors are limited. 

3.5.6 Scintillation 

Scintillation detectors, e.g., sodium iodine (Nal) or hquid xenon (LXe), are the solution 
to accumulate large mass of detector material 100 kg). However, it is more difficult 
to achieve radiopurity comparable to Ge detectors [51]. 

Moreover, as mentioned in Subsec. 3.5.2, scintillation detectors do not measure the 
elevation of the temperature in the crystal, but the light emitted by the electrons produced 
due to the ionization, thus the energy threshold for these detectors may be substantially 
higher than the thermal calorimeters [1]. 

Meanwhile, the Nal-bascd experiments, such as DAMA [48], NalAD [55], and EL- 
EGANT, originally attempted to use a pulse shape discrimination to statistically iden- 
tify the WIMP signals from their observed events (detailed discussions about the back- 
ground discrimination will be given in Sec. 3.6). It was found that the low number of 
detected scintillation photons per keV of incident energy (called "photo-electron per keV" 
or "p.e./keVee") restricts the usefulness of this method at low energy [51]. This means 
that the background may be problematic. The technique is now being investigated for Csl 
or CaF scintillator, where the difference in time constants between scintillations induced 
by electron- and nuclear- recoils are larger than in the Nal detectors [51]. 

3.5.7 Combinations of two different signals 

Actually, most of the direct WIMP detection experiments use detectors with mixed 
techniques and measure simultaneously two signals. For example, for cryogenic detec- 
tors, the CDMS and EDELWEISS collaborations [59], [60] investigate the heat-ionization 
signals (more details will be given in Subsecs. 3.7.1 and 3.7.4), and the CRESST col- 
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laboration [56] explores the heat-scintillation channel (more details will be given in Sub- 
sec. 3.7.2). 

Combining information measuring from two different channels can offer a powerful 
event-by-event rejection method for the background discrimination down to 5 to 10 keV 
recoil energy [37]. As mentioned in Subscc. 3.5.3, due to the quenching effect of the 
detector material, the ratio of the ionization or the scintillation signal to the heat signal is 
significantly different for the nuclear recoils and for the electron recoils. Similarly, nuclear 
recoils due to WIMP or neutron interactions have a much higher characteristic light 
over charge ratio than electron recoils due to electron and 7-ray interactions [61]. Thus 
simultaneous measurements of two of the heat, the ionization, or the scintillation signals 
can be used to distinguish nuclear recoils induced by WIMPs from electron recoils induced 
by electron or 7-ray interactions (more details about different methods for background 
discrimination will be given in Sec. 3.6). 

3.6 Background and background discrimination 

As mentioned in the beginning of this chapter, due to the very low cross section of 
WIMPs on ordinary material, the event number of the elastic WIMP-nucleus scattering 
is very rare and the backgrounds coming from different sources are much larger. 

For example, cosmic rays and cosmic-ray induced 7-rays with energies in the keV to 
MeV range, radioactive isotopes in and around the detector (in the equipment) should 
be considered. Moreover, neutrons induced by cosmic muons can produce nuclear recoil 
events similar to the real events induced by WIMPs. And electron recoils from photons 
(x-ray and 7-ray radiations) and electrons are also a major background. 

3.6.1 Cosmic muons and underground laboratories 

At ground level, approximately 10^ cosmic muons pass through per square centimeter 
of the Earth's surface per day [51]. They can induce nuclear transmutations to unstable 
isotopes throughout the detector volume [51]. 

In order to protect from the penetrating cosmic muon flux, it is necessary to place 
the detector in deep underground. In underground laboratories such as the Soudan 
Underground Laboratory (the CDMS collaboration) in Minnesota in the USA, the Gran 
Sasso National Laboratory (the CRESST and DAMA collaborations) in Italy, or the 
Laboratoire Souterrain de Modane (LSM, the EDELWEISS collaboration) in the Frejus 
Tunnel in the French-Italian Alps, the muon flux can be reduced by a factor of 10^ ~ 10^ 
[51]. 
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3.6.2 External natural radioactivity and (passive) shielding 

External sources of radioactivity mean the radioactive isotopes in the rock around the 
underground laboratory and in the walls of the laboratory. 

A shielding from external natural radioactivity can be achieved by surrounding the 
detector with thick absorbing material [51]: high-Z materials hke lead are very effective 
for stopping 7-rays with MeV energy, while \ow-Z materials are sufficient for stopping 
low energy 7-rays as well as a- and /3-radiations. 

3.6.3 Internal natural radioactivity and radiopure materials 

Beyond a thickness of 15 to 25 cm of lead shielding [51], one has to consider the 
internal radioactivity of the equipment, of the contamination near the detector or in the 
target material, and even of the lead shielding itself. 

Internal radioactivity can be reduced very well by using detectors (and the other 
experimental equipment) made of radiopure materials. Archeological lead has also often 
been used since it has already been shielded from cosmic rays for 2000 years. 

3.6.4 Active background rejection 

Except passive shielding around the detector, most experiments use also some different 
techniques for active background rejection. 

Generally there are two different types of background discrimination. Statistical re- 
jections are used to ascertain which fraction of the total event sample comes from a 
well-defined type of background, but cannot tell for one individual event [51]. Moreover, 
this kind of rejections depends strongly on the theoretical predictions about the true 
signals induced by WIMPs and the background events. 

On the other hand, the event-by-event rejections check each recorded event in the 
detector independently of the others ("blind") and can be used to reject background 
events with an almost 100% certainty. Note that, however, in practice there is always a 
small probabihty that some background event may fake the signals induced by WIMPs. 

3.6.5 Neutron induced nuclear recoils 

Cosmic muons can induce neutrons in the inner lead shielding and such fast neutrons 
can induce nuclear recoils similar to those induced by WIMPs. 

Fast neutron shielding consists of moderators made of material with a high density of 
hydrogen, such as polyethylene or water [51]. 

3.6.6 Multiple-scatter events and array of detectors 

The interaction between WIMPs and ordinary material is too weak, or, equivalently, 
the mean free path of a WIMP in ordinary matter is too long (of the order of a light-year 
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[51]), so that WIMPs could never interact more than once in a single detector or two 
adjacent detectors. In contrast, the mean free path of a neutron or a high energy photon 
is of the order of cm, thus multiple-scatter events produced by neutrons are more common 
[51]. 

Hence, an array of closely packed detectors (e.g., the tower with six detectors used by 
the CDMS experiment, see Subsec. 3.7.1 and their web page [59]) can efficiently identify 
these multiple-scatter events [51]. 

3.6.7 Electron recoils 

Theoretically WIMPs interact only with the nuclei (through the couphng to quarks) 
and produce nuclear recoils, while, due to the electromagnetic interaction, the dominant 
radioactive backgrounds interact usually with the electrons and produce electron recoils. 
Therefore, the experiments which can discriminate between the events due to nuclear 
recoils and events due to electron recoils can reject most radioactive background [62]. 

There are three ways to distinguish nuclear recoils from electron recoils. First, as 
mentioned in Subsec. 3.5.3, due to the quenching effect of the detector material, the ratio 
of the ionization or the scintillation signal to the heat signal is significantly different for 
nuclear recoils and for electron recoils. Thus one can measure simultaneously the heat 
signal and the ionization or the scintillation signal to distinguish the nuclear recoil events 
from the electron recoil events. 

Second, the decay times of pulses for nuclear recoils may be different than that for 
electron recoils [62] . Thus some experiments use only a scintillation detector, but measure 
also the timing of the signals [62] . However, due to the limited resolution and discrimina- 
tion power of this pulse shape analysis at low energies, this effect allows only a statistical 
background rejection [37]. This technique has been used by e.g., the DAMA/Nal and 
NalAD experiments (Nal(Tl) detector) [48], [55] and ZEPLIN-I experiment (hquid xenon 
detector) [55]. 

Third, as mentioned in Subsec. 3.5.3, the range of a electron recoil is of the order 
of lira and that of a nuclear recoil is only of the order of nm. Thus the nuclear recoils 
have a much larger energy loss per unit length, dQ/dx, or, equivalently, produce a much 
higher energy density, than the electron recoils. Therefore, some experiments are actually 
immune to electron recoils because the energy they deposit is not dense enough to trigger 
[62]. 

3.6.8 Surface events and self shielding 

Due to their very long mean free path, WIMPs will interact uniformly throughout the 
detector volume. In contrast, due to the short mean free path of the high-energy photons 
and neutrons, for a detector with a large volume, the interactions induced by radiations 
originating from the surrounding material and surface contamination will occur mostly 
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at the detector surface [51]. 

This "self-shielding" effect leads to the incentive of building large position-sensitive 
detectors in order to reject the surface events. Moreover, low-energy photons, a- and 
/3-rays have very short mean free path (< mm), and can be rejected even if the position 
sensitivity is limited [51]. 

3.6.9 Incomplete charge collection 

Some electromagnetic events occurring very near the detector surface can also mimic 
nuclear recoils because they produce less ionizations than expected from electron recoils 
[45]. But such surface events can also be rejected by the self-shielding effect. 

3.6.10 Shape of the recoil energy spectrum 

The shape of the recoil spectrum dR/dQ for a given WIMP mass in some simple halo 
models can be predicted numerically or even analytically (e.g., {dR/dQ)Gan in Eq.(3.21) 
and {dR/dQ)sh in Eq.(3.31)). The measured recoil spectrum should be consistent with 
the expectation. 

However, first, the overall shape of the expected spectrum is (approximately) expo- 
nential, as is the case for many background sources; second, different velocity distri- 
bution functions in different halo models could predict totally different recoil spectrum 
(cf. {dR/dQ)Ga.u in Eq.(3.21) and {dR/dQ)sh in Eq.(3.31)). Moreover, the expected signal 
events measured by the currently operated detectors and even next-generation ones are 
at most only a few per day. Hence, as we will see in Sec. 4.2, at present a meaningful 
reconstruction of the recoil spectrum with a small statistical error is actually impossible. 

3.7 Cyrogenic detectors 

As discussed above, a WIMP detector is constrained by three important requirements: 
low threshold, (ultra) low background, and high detector material mass [22]. 

In the following I will present some important collaborations worldwide and summarize 
their recent results and plans in the near future. More details about these collaborations 
and their experiments can be found in the references. 

3.7.1 CDMS 

The Cryogenic Dark Matter Search (CDMS) collaboration [59] uses the Berkeley 
Center for Particle Astrophysics (CfPA) germanium cryogenic detector [1]. Their first 
test run was at the Stanford University Underground Facility [1], and now moved to the 
deep Soudan Underground Laboratory (Soudan mine) in Minnesota in the USA [63]. The 
Soudan mine has 780 m rock overburden (2090 meters water equivalent, m.w.e.) [63], the 
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surface muon flux is then reduced by a factor of 5 x 10^ [63], and the neutron background 
is also reduced by a factor of 400 [22] (~ 4 x lO^^/kg/day [64]). 

It was the flrst experiment to operate a detector measuring simultaneously ionization 
and heat signals with a germanium/silicon crystal as the target material [22]. They 
developed Z(dcpth)-scnsitivc Ionization and Phonon (ZIP) detectors. The principle of 
their ZIP detector is basically the same as that discussed in Subsec. 3.5.5, except that 
the heat sensor is replaced by a thin film sensor and thus able to detect phonons before 
their complete thermalization [51]. 

Their "tower (s)" with mixed Ge and Si detectors are powerful for subtracting the 
neutron background [63]. Except the neutron multiple-scatter events discussed in Sub- 
sec. 3.6.6, while Ge and Si have similar scattering rates per nucleon for neutrons, the 
WIMP-nucleon scattering rate is expected to be 5-7 times greater in Ge than in Si for 
all but the lowest-mass WIMPs. Moreover, the kinematics of neutron elastic scattering 
gives a recoil energy spectrum scaled in energy by a factor of ~ 2 in Si compared to 
Ge, whereas the factor would be ~ 1 or less for WIMP elastic scattering. All of these 
three methods can be used (together), in conjunction with Monte Carlo simulations, to 
statistically subtract any neutron background. 

In addition, because the athermal phonons from electron recoils are faster than those 
from nuclear recoils, particularly if the electron recoils occur near a detector surface, 
by collecting such fast, athermal phonons with thousands of thin-film sensors, their ZIP 
detector can discriminate very well against the surface electron recoils [62]. 

According to Ref. [63], for recoil energies above 10 keV, events due to the background 
photons can be almost perfectly rejected, and more than 96% of the incomplete charge 
collection events can also be rejected by using additional information from the shape, 
timing, and energy partition of the phonon pulses (namely, only events with both slow 
phonon pulses and low ionization have been accepted) , while over half of the nuclear- recoil 
events should be kept [62]. 

In the first Soudan run of the CDMS-II experiment (from October 11, 2003 to January 
11, 2004) [63], one tower with 4 Ge (each 250 g) and 2 Si (each 100 g) ZIP detectors has 
been operated for 52.6 live days, the recoil energy thresholds of these six detectors were 
between 10 and 20 keV, only one candidate event with a recoil energy of 64 keV in one 
Ge detector has been measured. For a WIMP mass of 60 GeV/c^, a 4 x 10~^ pb upper 
limit on the spin-independent WIMP-nucleon cross section ^ from Ge as been achieved. 
Meanwhile, thanks to the ^^Ge (spin-|) and ^^Si (spin-i) content of natural germanium 
and silicon, a 0.2 pb upper hmit on the spin-dependent WIMP-neutron cross section 
for a WIMP mass of 50 GeV/c^ has also been achieved. These were the world's lowest 
limits on the WIMP-nucleon cross-section in the case of spin-independent interactions 
and spin-dependent interactions with neutrons. 

^Here the cross-section ctq shown in Eqs.(3.5) and (3.13) is normaUzed to a single nucleon a^n in order 
to allow comparisons between different target nuclei. 



56 



In the second Soudan run (from March 25 to August 8, 2004) [65], [66], and [46], 
two towers (one tower with 4 Ge and 2 Si detectors and the other one with 2 Ge and 
4 Si detectors) have been operated for 74.5 hve days and the recoil energy thresholds 
have been improved to be only 7 keV, one more candidate event with a recoil energy of 
10.5 keV in one Ge detector has been measured. For a WIMP mass of 60 GeV/c^, the 
upper limit on the spin-indcpcndcnt WIMP-nucleon cross section has been given as 1.6 
xlO~^ pb from Ge and 3.4 xlO^® pb from Si (see Fig. 3.4). These limits are a factor 
of 6 lower than those given by the ZEPLlN-1 experiment [55] (sec Subscc. 3.8.5) and an 
order of magnitude lower than those of the CRESST and EDELWEISS collaborations 
[62]. Moreover, their results excluded the overlap between the CDMS and DAMA/Nal 
allowed regions at WIMP masses > 25 GeV/c^, though compatible regions at lower masses 
remain [66] (see Fig. 3.4). 

Now the CDMS collaboration is preparing for five towers with totally 19 Ge (4.75 kg) 
and 11 Si (1.1 kg) ZIP detectors, and will improve their sensitivity a factor of ^ 10 [46]. 
Furthermore, they also planned a SuperCDMS project which will start with a total mass 
of 25 kg (Phase A, each detector will be 640 g) and be improved to 150 kg (Phase B) and 
eventually 1000 kg (Phase C) [46] , in order to achieve ~ 10~^ pb sensitivity (for a WIMP 
mass of 60 GeV/c^, Phase A, see Fig. 3.5) [46], corresponding to C(10~^) events/kg/day 
event rate in the energy range between 15 and 45 keV [64]. They will also move to 
the SNO Underground Laboratory at the Sudbury mine in Canada. The ~ 6000 m.w.e. 
overburden at this site results in over two orders of magnitude suppression in the neutron 
background compared to Soudan [64]. 

3.7.2 CRESST 

The Cryogenic Rare Event Search using Superconducting Phase Transition Thermome- 
ters (CRESST) collaboration [56] uses heat-scintillation detectors with CaW04 crystal 
in the Gran Sasso National Laboratory in Italy. Their detector provides a good rejection 
of surface events as of photons due to the much larger light yield from all electron recoils 
relative to nuclear recoils [62] . 

As mentioned in Subsec. 3.5.4, their detector uses the superconducting-normal phase 
transition due to the difference of the temperature. A thin superconducting film of tung- 
sten (W) has be grown on a silicon detector and held just below the critical temperature. 
Heat produced by WIMP-nucleus scatterings will change the film to its normal state and 
the change in resistance could be measured [1]. However, as mentioned in Subscc. 3.5.6, 
the threshold energy of such a scintillation detector is relatively higher than for an ion- 
ization detector. Thus a disadvantage of the CRESST heat-scintillation detector is that 
an event measured by the phonon channel but producing no light may mimic a WIMP 
signal [62]. 

In 2003 CRESST ran two prototype detectors for a couple of months without neutron 
shielding. A significant neutron background on the oxygen in their CaW04 detectors was 
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observed [62] and the light yield for W recoils is significantly less than for Ca or O recoils. 
This result indicates that WIMPs are expected to interact primarily with W nuclei, while 
neutrons will interact relatively more often with O and Ca nuclei [51]. 

In early 2004 they operated two 300 g CRESST-II prototype detector modules, 16 
events have been recorded and a rate for nuclear recoil energy between 12 and 40 keV of 
(0.87±0.22) events/kg/day has been obtained [67]. However, this is compatible with the 
rate expected from neutron background, and most of these events lie in the region of the 
phonon-light plane anticipated for neutron-induced recoils [67]. Moreover, a particularly 
strong limit for WIMPs with coherent scattering results from selecting a region of the 
phonon-light plane corresponding to tungsten recoils, where the best module shows zero 
events [67]. The sensitivity achieved by the CRESST-II experiment is given in Fig. 3.4. 

Now they are preparing the scientific run of the CRESST-II experiment with 33 
detectors (each 300 g) and totally ~ 10 kg target material (see Fig. 3.5). 

3.7.3 DAMA 

The DArk MAtter (DAMA) collaboration [48] uses a scintillation detector with ~ 100 
kg Nal(Tl) in the Gran Sasso National Laboratory (Laboratori Nazionali del Gran Sasso, 
LNGS) in Italy [52]. With 1400 m rock overburden (3500 m.w.e.), the total muon fiux is 
reduced to ~ 1/m^/hr (one order of magnitude lower than that of CDMS), the external 
7-ray flux is reduced by a factor of 10^. 

They are the only collaboration which claimed to detect the signal of halo Dark 
Matter due to the annual modulation effect discussed in Sec. 3.2. Figs. 3.3 show the 4- 
year and the 7-year results of the DAMA/Nal experiment [53], [54], their threshold energy 
is about 2 keVge, corresponding to approximately 22 keV recoil energy [51]. Meanwhile, 
they pubhshed a WIMP mass m-^ ~ 52 GeV/c^ and a WIMP-proton cross section a^,^ ~ 
7.2 X 10~^ pb [52] under the standard assumptions of WIMP halo described in Subsec. 1.3.1 
(see Fig. 3.4). 

However, the DAMA collaboration uses the pulse shape discrimination (PSD) tech- 
nique (see Subsec. 3.6.7) to statistically (not event-by-event) discriminate the measured 
events [51]. On the other hand, the CDMS results [68], [63], [65] and [66] are clearly 
incompatible with the signal claimed by DAMA under the standard assumptions of the 
WIMP halo and spin-independent WIMP-nucleus coupling [62]. 

But these two experiments might still be compatible in some exotic scenarios. One 
possibility is to postulate rather light (m^ < 10 GeV/c^) and fast WIMPs with large 
scattering cross section [69]. Since for this case mj ^ m^, the maximal recoil energy 
induced by the scattering on iodine will then be smaller than the threshold energy of the 
DAMA's detector and the recorded events have thus been induced by the scattering on 
sodium (see Eq.(3.46)). Similarly, the Ge nuclei used by the CDMS experiment could 
also be too heavy to deposit recoil energy large enough to be measured. However, the 
null results from the Si detector used by CDMS (see Subsec. 3.7.1) should (almost) rule 
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Figure 3.3: The published resuhs of the DAMA/Nal experiment, (a) In the 2-6 keVge 
cumulative energy interval over 4 annual cycles, since January 1st of the first year of data 
taking. Theoretically expected minimum (dashed line), maximum (dotted line) (figure 
from [53]). (b) In the 2-4 keVge, 2-5 keVee, and 2-6 keVee cumulative energy intervals 
over 7 annual cycles and end of data taking in July 2002 (figure from [54] ) . 
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out this possibility. Another possible way out is to postulate that the detected events are 
actually inelastic, leading to the production of a second particle that is almost, but not 
exactly, degenerate with the WIMP [70] . 

They are running now the DAMA/LIBRA (Large sodium Iodide Bulk for RAre pro- 
cesses) experiment with totally ~ 250 kg Nal(Tl) [53]. 

3.7.4 EDELWEISS 

The Experience pour DEtecter Les WIMPs En Site Souterrain (EDELWEISS, EDW) 
collaboration [60] is in the Laboratoire Souterrain de Modane (LSM) in the Prejus Tunnel 
in the Prench-Italian Alps. With ~ 1800 m rock overburden (~ 4800 m.w.e.), the muon 
flux can be reduced to ~ 4/m^/day, the fast neutron flux can be reduced to ~ 1.6 x 
10-VcmVs [71]. 

The EDELWEISS-I experiment has used 3 (each 320 g) cryogenic heat-and-ionization 
Ge detector [51]. Their calibrations indicate that the larger ionization/recoil energy ratio 
of electron recoils results in very good discrimination ability against photon backgrounds 
down to the 20 keV threshold energy [62]. 

In 2000 to 2002 [47], the EDELWEISS-I experiment has been operated for an exposure 
of 13.6 kg-day. The energy threshold was 13 kcV and no event has been recorded. In 2003 
[47], the second run of the EDELWEISS-I has been operated for an exposure of 48.4 kg- 
day (totally 62 kg-day) and 40 nuclear recoil candidate events have been recorded in the 
energy range 15 to 200 keV: 18 events between 15 and 20 keV, 16 events between 20 and 
30 keV, 3 events between 30 and 100 keV, 3 events between 100 and 200 keV, and more 
19 events have been observed below 15 keV; most likely due to remaining background 
neutrons and surface electrons. According to these results, they gave exclusion limits for 
WIMP masses above 25 GeV/c^ [47] (see Fig. 3.4). 

Now they are preparing the EDELWEISS-II experiment with 120 detectors [71] (see 
Fig. 3.5). 

3.7.5 Heidelberg-Moscow (HDMS) 

The Heidelberg-Moscow collaboration [58] uses a ~ 2 kg ''^Ge semiconductor ioniza- 
tion detector in the Gran Sasso National Laboratory in Italy and achieved a very low 
background count rate (< 0.2 event/kg/day) in the interval 10 ~ 40 keV, and a threshold 
energy Qthre 4 ~ 10 keVge (equivalent to ~ 15 ~ 30 keV recoil energy) [22]. 

They produced the flrst limits on WIMP searches and until recently had the best 
performance. However, without positive identiflcation of nuclear recoil events, these 
experiments could only set limits, e.g., excluding sneutrinos as major component of the 
galactic halo [37]. 
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3.7.6 KIMS 



The Korea Invisible Mass Search (KIMS) collaboration [72] in the Yangyang Under- 
ground Laboratory (Y2L, ~ 700 m rock overburden) in South Korea developed a CsI(Tl) 
crystal scintillation detector and uses an improved pulse shape discrimination to discrimi- 
nate nuclear recoil events induced by WIMPs from electron recoil events induced by 7-ray 
background [73]. 

They have operated 4 detectors with an exposure of 3407 kg-day in the temperature 0° 
C and null signals have been observed [73]. Due to their ^^^Cs and ^^^I target nuclei this 
result has been used to give the lowest upper hmits on the spin-dependent WIMP-proton 
cross section for WIMP masses > 30 GeV/c^ [73] (the lowest upper hmit on the spin- 
dependent WIMP-neutron cross section has been obtained by the CDMS collaboration, 
see Subsec. 3.7.1). Moreover, although several experiments have already given exclusion 
limits rejecting the DAMA signal region, it is the first time that such a limit obtained by 
the crystal scintillation detector containing ^^^I, which has been usually assumed to be 
the dominant nucleus for the spin-independent interaction in the Nal(Tl) crystal. 

3.7.7 PICO-LON 

The Planar Inorganic Crystals Observatory for LOw-background Neutr(al)ino (PICO- 
LON) collaboration in Japan uses detector with multilayer (3 layers for PICO-LON-0 
and 16 layers for PICO-LON-I) thin Nal(Tl) crystals [74], [75]. 

A special advantage of their 0.05 cm thin and 5 cm x 5 cm wide area Nal(Tl) crystals 
is the position sensitivity of the recoil events [75]. The position resolution for the thinner 
directions is as good as 0.05 cm due to the segmentation of the detector [75]; while, the 
position information in the wider area was obtained by analyzing the ratio of the number 
of photons collected on the opposite sides of the detector [75] and a circle with ~ 0.5 cm 
radius has been reached [76]. Moreover, they have also claimed that a very low threshold 
energy ~ 2 kcVge has been measured [74], [76]. 

The PICO-LON-0 experiment has been run at the surface laboratory at Tokushima, 
and the PICO-LON detector will be installed into Oto Cosmo Observatory (100 km south 
from Osaka) covered by thick rock with ~ 1200 m.w.e. [75]. 

3.8 Liquid noble gas detectors 

Liquid noble gas detectors have the advantage of easier scaling to large masses since 
it is based on hquids [62]. They can also be allowed to operate in higher temperatures: 
165 K for xenon, 87 K for argon, and 27 K for neon [77]. 

Due to its high-A value, liquid xenon (LXe) has been used by the ZEPLIN collabora- 
tion [55] as the first liquid noble gas detector material (more details about the ZEPLIN 
experiments will be given Subsec. 3.8.5). Recoils in the liquid noble gas such as Xe can 



61 



induce both ionization and excitation of Xe atoms. An excitation produced by a nuclear 
recoil usually induces emission with a single photon, whereas that reduced by an electron 
recoil emits photons in form of a slow triplet, thus nuclear recoils have a faster pulse 
shape than electron recoils [62] , [78] . 

Besides xenon, argon and neon are also suitable detector materials. The effect of 
faster pulse shape is particularly clear in Ar and Ne, leading to their extremely good 
discrimination based on timing [62]. In addition, due to the form factor effect introduced 
in Eq.(3.5) (two nuclear form factors for spin-independent cross section have been given 
in Subsec. 3.1.2), the event rate in e.g., argon is less sensitive to the energy threshold 
than in xenon [61]. 

Moreover, as discussed in Subsec. 3.5.7, the scintillation (light) over ionization (charge) 
ratio can be used additionally to discriminate the nuclear recoils from the electron recoils 
due to electron and 7-ray interactions. 

However, discrimination of the radioactive contamination in the detector material, 
such as ^^Kr in hquid Xe (25 ppm Kr in natural Xe) or especially ^^Ar in liquid Ar [62] 
as well as of the surface radioactivity from the liquid container [51], and relatively larger 
threshold energies could be primary challenges for detectors using liquid noble gases. 

3.8.1 ArDM 

The Argon Dark Matter (ArDM) experiment at CERN [79] uses a ton-scale detec- 
tor with liquid argon (LAr), which measures simultaneously the scintillation and the 
ionization signals [80], [61], and [81]. 

With a recoil energy threshold of 30 keV and a sensitivity of 10~^ pb WIMP-nucleon 
cross section, the ArDM experiment has been expected to yield approximately 100 events 
per day per ton [61]. By improving the background rejection power and further limiting 
the background sources, a sensitivity of 10~^ pb (1 event per day per ton) would become 
reachable [61]. 

3.8.2 WARP 

Similar to the ArDM experiment, the WIMP ARgon Programme (WARP) experiment 
[82] uses also a dual-phase (gas-liquid) argon detector [62] in the Gran Sasso Underground 
Laboratory in Italy [77] . By using a strong electric field, ionization electrons are drifted 
out of the liquid argon into gaseous phase, where they are detected via the secondary 
luminescence [51], [62]. The discrimination technique is based on both of the pulse shape 
of the photon emissions and on the ratio of the scintillation to ionization energies [62] 
described in Subsec. 3.5.7. 

Their first run of a 2.6 kg (1.87 t) prototype with a 96.5 kg-day exposure resulted in 
no candidate events above the threshold energy 55 keV [77] (see Fig. 3.4). Later they will 
upgrade the detector to totally 140 kg (> 100 £) [77]. 
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3.8.3 XENON 



The XENON collaboration [83] uses also a dual-phase xenon time projection chamber 
(XeTPC) with 3D position sensitivity in the Gran Sasso Underground Laboratory in Italy 
[84]. 

The XENONIO experiment uses 15 kg liquid xenon, the background rate from ^^Kr 
contamination is reduced by a factor of 5000 by using a commercially available low-Kr (5 
ppb) xenon, the self-shielding effect (see Subsec. 3.6.8) of LXe can reduce the background 
events in the central region of the LXe target by more than one order of magnitude for 
recoil energy below 50 keVee, (5 to 15 keVee corresponds to roughly 15 to 40 keV recoil 
energy) [84]. 

By measuring simultaneously the scintillation and the ionization produced by radia- 
tion in pure liquid xenon, the detector can discriminate signals from background down to 
4.5 keVee [85]. Between October 6, 2006 and February 14, 2007 the XENON experiment 
has been run for 58.6 live days and 10 events have been observed in the energy range 4.5 
to 26.9 keVee- However, none of these events are likely WIMP interactions. [85] 

Their newest result gives a 90% C. L. upper limit on the spin-independent WIMP- 
nucleon cross section of 8.8 x 10~^ pb for a WIMP mass of 100 GeV/c^ [85], a factor 
of 2.3 lower than the limit achieved by CDMS-II experiment (see Subsec. 3.7.1). For a 
WIMP mass of 30 GeV/c^, the limit is 4.5 x 10"^ pb [85] (see Fig. 3.4). 

The XENONIO experiment will be upgraded to ~ 100 kg [51] and a WIMP-nucleon 
sensitivity of 2 x 10"^ pb in 2008 [86] (see Fig. 3.5). 

3.8.4 XMASS 

The Xenon Neutrino MASS Detector (XMASS) experiment uses a 100 kg (intend 
to ultimately 800 kg) [37] single-phase Xe detector [62] at the SuperKamiokande site in 
Japan. 

Because of their 100 kg total mass of target material, XMASS has a good position 
sensitivity and has demonstrated the self-shielding effect (see Subsec. 3.6.8) to reduce 
the background events induced by multiple scattering (see Subsec. 3.6.6) and surface 
contamination [37]. 

Besides the XMASS experiment with hquid Xe, CLEAN (Cryogenic Low-Energy As- 
trophysics with Neon) and DEAP are also single-phase experiments. They use Ar or 
Ne as detector material in oder to take advantage of the much larger timing difference 
between nuclear recoils and electron recoils described above [62], [87]. 

3.8.5 ZEPLIN 

The Zoned Proportional Scintillation in Liquid Noble Gases (ZEPLIN) collaboration 
[55] first used a liquid xenon scintillation detector in the Boulby Laboratory (1070 m 
underground) in the United Kingdom. 
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Figure 3.4: The curves show the sensitivities of (the exclusion limits on) the spin- 
independent WlMP-nucleon cross section versus WIMP mass achieved by the CDMS 
[65] (the blue solid and the blue dashed lines, for Ge and Si, respectively), the CRESST 
[67] (the black solid line), the DAMA [52] (the green area), the EDELWEISS [47] (the 
red solid line), the WARP [77] (the dark green solid line), the XENON [85] (the black 
dashed line), and the ZEPLIN [78] (the red dashed line) collaborations (plot generated 
by http : // dmtools . berkeley . edu/limitplots/). 

ZEPLlN-1 was a single-phase experiment with 6 kg (3.1 kg fiducial) Xe, so it could 
not collect the ionization signals, but depended solely on the pulse shape discrimination 
[62]. In an exposure of 293 kg-day, no excess consistent with nuclear recoils was seen 
[62]. However, the published limits are somewhat controversial, because their calibration 
results of the neutron recoil discrimination do not appear to be convincing enough to 
consider the limits set on the WIMP signal to be as reliable as the ones set by the 
cryogenic experiments [37]. Actually, with only 1.5 photo-electron per keV and a three- 
fold coincidence, searching for the WIMP signal in the 2-10 keVee region is for ZEPLlN-1 
quite challenging [37]. 

ZEPLlN-11 has been upgraded to a dual-phase experiment with 31 kg Xe [78]. In the 
first run with 225 kg-day exposure, 29 events have been observed in an acceptance window 
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defined between 5 keVee and 20 keVge- With a summed expectation of 28.6 ± 4.3 7-ray 
and radon progeny induced background events, these figures provide a 90% C. L. upper 
hmit to the number of nuclear recoils of 10. 4 events in this acceptance window, which 
converts to a spin-independent WIMP-nucleon cross-section of 6.6 x 10"^ pb for a WIMP 
mass of 65 GeV/c^ (see Fig. 3.4). For the second run a sensitivity of 2 x 10~^ pb has 
been expected (see Fig. 3.5). 

3.9 Superheated droplet detectors (SDD) 

Metastable liquid droplets immersed in a gel expand (explode) due to a phase transition 
to the gaseous phase, when a particle or nucleus with sufficiently high energy density 
(energy deposited over unit length, dQ/dx) interacts in the liquid [45]. A broad range of 
detector materials could be used and inexpensively scaled to large masses [62]. 

The main (best) advantage of such integrating detectors is that, by tuning pressure 
and temperature, the threshold energy of the detector can be adjusted so that the detector 
could be insensitive to the low energy density [62] . Thus electron recoils and a-radiation 
events can be rejected automatically. 

3.9.1 DRIFT 

The Directional Recoil Identification From Tracks (DRIFT) experiment [55] uses a low 
pressure Xe-CS2 gas mixture TPC [22]. Using the negative CS2 ions instead of usual e~ 
as charge carriers can reduce the diffusion and thus achieve millimetric track resolution 
[22] . The ionization tracks will be measured with a multi-wire proportional chamber in a 
low-pressure gas [45], this offers the most convincing technique to measure the direction 
of nuclear recoils [37], [87]. 

However, one disadvantage of the DRIFT'S detector is the very low target mass and/or 
the need of a huge detector volume. 

3.9.2 MIMAC-He3 

The MIcro-tpc (temporal projection chambers) MAtrix of Chambers of Helium 3 
(MIMAC-He3) experiment [88] uses an ultra cold pure ^He detector [37]. The use of ^He 
as target material is motivated by its privileged features for Dark Matter search compared 
with other target nuclei. First, ^He being a spin-| nucleus with a single neutron, a detector 
made of such material will be sensitive to the "neutron spin-dependent" interaction, 
leading to a natural complementarity to most existing or planned Dark Matter detectors 
as well as proton based spin-dependant detectors [90]. Moreover, the mass of the ^He 
atom is 2.81 GeV/c^, the recoil energy range is expected less than 10 keV [89]. Hence, 
MIMAC-He3 could be used to measure events with low recoil energies as well as detect 
low-mass WIMPs (6 GeV/c^ < m^^ < 40 GeV/c^ [90]). 
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There are several more advantages for a detector using ^He [89], [90]: first, there are 
no intrinsic x-rays; second, a very low Compton cross section to 7-ray (two orders of 
magnitude weaker than in Ge) can reduce the natural radioactive background by several 
orders of magnitude; third, the neutron signature due to the capture process: n + ^He 
— ^ p + ''H + 764 keV, will be very useful for discrimination of neutron background. On 
the other hand, the double detection of the ionization energy and the track projection by 
means of the TPC chambers can assure the electron-recoil discrimination. 

However, similar to the DRIFT experiment, one disadvantage of their detector is the 
very low target mass and/or the need of a huge detector volume. 

3.9.3 PICASSO 

The Project In CAnada to Search for Supersymmetric Objects (PICASSO) experiment 
in the Sudbury Neutron Observatory (SNO) in Canada uses ■'^^F (spin-| isotope) as detec- 
tor material and search for spin-dependent Dark Matter particles (see Subscc. 3.4.2) [22]. 

The principle of the PICASSO detector is as follows: Small superheated freon droplets 
imbedded in a gel matrix at room temperature. The nuclear recoil of ^^F induces the 
explosion of a droplet, creating an acoustic shock wave which will be measured with 
piezoelectric transducers [22]. 

3.9.4 SIMPLE 

The Superheated Instrument for Massive ParticLe Experiments (SIMPLE) collabo- 
ration in the Laboratoire Souterrain a Bas Bruit (LSBB, ~ 1500 m.w.e.) in Frence uses 
C2CIF5 and CF3I and searches also for spin-dependent Dark Matter interaction [91], [92]. 
Their results exclude the spin-dependent WIMP-proton cross section above 1.14 pb for a 
WIMP mass of 50 GeV/c^ [93] (an upper limit of 0.2 pb on the spin-dependent WIMP- 
neutron cross section has been given by the CDMS collaboration, see Subsec. 3.7.1). 

3.10 Prospects 

So far we did not obtain any convincing signal from experiments searching for Dark 
Matter particles. In the future, detector technique, better sensitivities as well as better 
background discrimination, will be improved. We need also some new ideas for detector 
building as well as application of experimental data. 

As described in Subsec. 3.7.1, the CDMS collaboration has achieved a (so far the best) 
sensitivity of ~ 10~^ pb for spin-independent WIMP-nucleon cross section and of ~ 10~^ 
pb for spin-dependent WIMP-neutron cross section. Together with the other collabora- 
tions described above, direct WIMP detection experiments have started to probe some 
possible regions in the parameter space predicted by some supersymmetric models. For 
next-generation detectors, sensitivities will be upgraded down to ~ 10~^ pb (see Fig. 3.5) 
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Figure 3.5: The curves show the sensitivities of the spin-independent WIMP-nucleon 
cross section versus WIMP mass projected by the SuperCDMS (25 kg), the CRESST- 
II (with CaW04), the EDELWEISS-II, the XENONIOO (100 kg), and the ZEPLIN- 
II experiments. The indications of the lines as in Fig. 3.4 (plot generated by 
http : / /dmtools . berkeley . edu/limitplots/). 

and, in long term, even ~ 10^^° pb, and the corresponding WIMP-nucleus scattering 
event rate is then ~ 5 events/ton/yr for Ge [37], as needed to probe large regions of 
MSSM parameter space. The total mass of detector material will also be improved to the 
ton scale. For example, the ~ 100 kg Xe detector of the XENON and XMASS collabo- 
rations. The CDMS collaboration is also preparing for their SuperCDMS projects with 
maximum 1100 kg target mass, while the CRESST and the EDELWEISS collaborations 
will also build to a larger collaboration "EURECA" (European Underground Rare Event 
search with Calorimeter Array). 

A nearly perfect background discrimination capability for next-generation detectors is 
also necessary. The ultimate neutron background will only be identified by the multiple 
interactions in a finely segmented or multiple interaction sensitive detector, and/or by 
operating detectors containing different target materials within the same setup [37]. 
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Furthermore, the measurement of the recoil directions of the target nuclei would 
provide additional information on the distribution of WIMPs in our Galaxy [62]. As 
described in Subsec. 3.9.1, the most convincing way for determining the recoil direction 
is by drifting negative ions in a temporal projection chamber to accurately record the 
tiny recoil distance. The DRIFT experiment has provided a proof of the principle, but 
it remains to be seen if such gas detectors with enough target material can detect some 
signals [62]. 

By the way, in order to present the WIMP-nucleon cross sections and the detector 
sensitivities in the future more conveniently and also suitably, we may consider to use 
"zepto" (10^21) or even "yocto" (IQ-^^) [94] barn instead of lO'^ pb (10''^^ cm^) or IQ-^^ 
pb cm2). 
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Chapter 4 



Reconstruction of the Velocity 
Distribution of WIMPs 



So far most theoretical analyses of direct detection of halo WIMPs, as discussed in 
Subsecs. 3.1.3 and 3.2.1, have predicted the detection rate for a given (class of) WIMPs, 
based on a specific model of the galactic halo (e.g., [24], [26], and [95]). The goal of my 
work is to invert this process. That is, I wish to study, as model-independently as possible, 
what future direct detection experiments can teach us about the WIMP halo. 

In this chapter I will use a time-averaged recoil spectrum, and assume that no di- 
rectional information exists. One can thus only reconstruct the (time-averaged) one- 
dimensional velocity distribution, fi(v). In the first section I will show how to derive 
the velocity distribution of WIMPs from the functional form of the recoil spectrum; the 
assumption here is that this functional form has been determined by fitting the data 
of some (future) experiment (s). I will also derive formulae for moments of the velocity 
distribution function, such as the mean velocity and the velocity dispersion of WIMPs, 
which can be compared with model predictions. 

Then I will present the method for reconstructing the velocity distribution of WIMPs 
directly from recorded signal events. This allows statistically meaningful tests of predicted 
distribution functions. I will also discuss how to estimate the moments of the velocity 
distribution directly from these data. 

Finally, I will show how to determine the mass of halo WIMPs, which one needs 
for the reconstruction of (the moments of) the velocity distribution, from experimental 
data directly. This allows also a useful comparison of the detected WIMPs with the new 
particle(s) produced at coUiders, e.g., the Large Hadron CoUider (LHC). 
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4.1 From the scattering spectrum 



In this section I will start with the differential rate for elastic WIMP-nuclues scattering 
given in Eq.(3.12) 

dR 



dQ 



AF'iQ) 



poo 






V 



dv , 



(3.12) 



and show how to find an expression for the one-dimensional velocity distribution function 
fi{v) for given (as yet only hypothetical) measured recoil spectrum dR/dQ. To that end, 
I first define 

dF,{v) Mv) 



dv V ' 

i.e., Fi{v) is the primitive of fi{v)/v. Then Eq.(3.12) can be rewritten as 
1 
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(4.1) 
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Since WIMPs (as candidate for CDM) in today's Universe move quite slowly, fi{v) must 
vanish as v approaches infinity. 



h{v^oo)^0. 



Hence 
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This means that Fi{v — > oo) approaches a finite value. ^ Differentiating both sides of 
Eq.(4.2) with respect to Vmm and using Eq.(3.10), one can find that (detailed derivations 
will be given in App. A. 2) 
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Since this expression holds for arbitrary I'min, one can write down the following result 
directly: 
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The right-hand side of Eq.(4.6) depends on the as yet unknown constant A. Recall, 
however, that fi{v) is the normalized yelocitj distribution, i.e., it is defined to satisfy 
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^The other properties of Fi{v) will be discussed in App. A.l. 
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Therefore, the normahzed one-dimensional velocity distribution function can be expressed 
as 



h{v)=J\f\-2Q 



d 

dQ 



1 fdR\ 



(4. 



with the normalization constant J\f (a detailed derivation will be given in App. A. 3) 
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Note that the integral here starts at Q = 0. 

In the next step I want to compute the moments of the velocity distribution function: 



1 (Qthre) 



v''fi{v)dv. 



(4.10) 



Here I have introduced a threshold energy Qthre- This is needed experimentally, since at 
very low recoil energies, the signal is swamped by electronic noise. Moreover, later we 
will meet expressions that (formally) diverge as Q — > 0. t'min(Qthre) is calculated as in 
Eq.(3.10). Inserting Eq.(4.8) into Eq.(4.10) and integrating by parts, one can find that 
(a detailed derivation will be given in App. A. 3) 
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Physically, {v) is the average WIMP velocity, while (w^) gives the velocity dispersion. ^ 
One emphasis here is that Eqs.(4-ll) and (4-12) can he evaluated directly once the recoil 
spectrum is known; knowledge of the functional form of fi{v) in Eq.(4.8) is not required. 
Note that the first term in Eq.(4.11) vanishes for n > if Qthre 0. In the same limit, 
(v^) J\faIo{0)/2 ^ 1 by virtue of Eq.(4.9). On the other hand, as written in Eqs.(4.8) 
and (4.9), the velocity distribution integrated over the experimentally accessible range 
of WIMP velocities gives a value smaller than unity. Using only quantities that can 
be measured in the presence of a non-vanishing energy threshold Qthrej Eq.(4.9) can be 
replaced by 



A/'(Qthre) = - 
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— -I- /o(Qthre) 
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(4.13) 



Using Af (Qthre) in Eq.(4.8) ensures that the velocity distribution integrated over v > 
^min (Qthre) gives unity. 

^The dispersion of the function fi{v) in the statistical sense is given by (w^) — (w)^. 
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Prom Eqs.(4.8), (4.9), and (4.11) to (4.13), it can be found that the (unreahstic) 
"reduced" spectrum (i.e., the recoil spectrum divided by the squared nuclear form factor) 
is more useful than the recoil spectrum itself. Meanwhile, note that the reduced spectrum 
from {dR/ dQ)G^.u in Eq.(3.21): 



is exactly exponential. This remains approximately true for the potentially quasi-realistic 
spectrum in Eq.(3.31) as well: 



In order to test these formulae, I have substituted the spectra in Eqs.(3.21') and (3.31') 
into Eqs.(4.8) and (4.11), taking Qthre = 0. They reproduced the normalized velocity 
distribution functions in Eqs.(3.20) and (3.29), as well as the results in Eqs.(3.22), (3.23), 
(3.32), and (3.33). The detailed calculations will be given in Apps. B.2.1 and B.2.2, 

respectively. 

One emphasis here is that the final results in Eqs.(4.8) and (4-11) independent 
of the as yet unknown quantity A defined in Eq.(3.13). They do, however, depend on 
the WIMP mass through the coefficient a defined in Eq.(3.11). This mass can be 
extracted from a single recoil spectrum only if one makes some assumptions about the 
velocity distribution fi{v). In the kind of model-independent analysis pursued here, 
can be determined by requiring that the recoil spectra using two different target nuclei 
lead to the same moments of the velocity distribution, (v'*), through Eq.(4.11). This 
method will be discussed in Sec. 4.3. Note that this can be done independent of the 
detailed particle physics model, which determine the value of ctq for two target nuclei. But, 
one will need to know both form factors of the target nuclei, which strongly depend on 
whether spin-dependent or spin-independent interactions dominate (see Subsecs. 3.4.1 to 
3.4.3). On the other hand, within a given particle physics model, the best determination 
of rriy^ should eventually come from experiments at high-energy colliders. However, we 
need also an alternative method as cross-check to prove whether the particle produced 
at coUiders is the same particle detected by the direct WIMP detection. 

4.2 From experimental data directly 

In the previous section I have derived formulae for the normalized one-dimensional 
velocity distribution function of WIMPs, fi{v), and for its moments (f"), from the recoil- 
energy spectrum, dR/dQ. In order to use these expressions, one would need a functional 
form for dR/dQ. In practice this might result from a fit to experimental data. Note that 
the expression for fi{v) in Eq.(4.8) requires knowledge not only of dR/dQ, but also of 




(3.21') 




(3.31') 
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its derivative with respect to Q, i.e., we need to know both the spectrum and its slope. 
This will complicate the error analysis considerably, if dR/dQ is the result of a fit. 

In this section I therefore go one step further, and derive expressions that allow 
to reconstruct fi{v) and its moments directly from the data. The assumption we have 
to make is that the sample to be analyzed only contains signal events, i.e., is free of 
background. This should be possible for the modern (next-generation) detectors (detailed 
discussions about the background and its discrimination have been given in Sec. 3.6). 
Having a sample of pure signal events, we can proceed with a complete statistical analysis 
of the precision with which we can reconstruct /i(f ) and its moments. 

However, in the absence of a true experimental sample of this kind, 1 had to resort to 
Monte Carlo experiments with an unweighted event generator written by M. Drees. Since 
detectors without directional sensitivity have been assumed, a single event is uniquely 
characterized by the measured recoil energy Q. Existing experiments such as CDMS 
[59] and CRESST [56] can determine the recoil energy quite accurately (some details 
about their detectors and experiments have been given in Subsecs. 3.7.1 and 3.7.2). We 
will see later that the statistical errors on the reconstructed velocity distribution fi{v) 
will be quite sizable even for next-generation experiments, given existing bounds on the 
scattering rate. It should therefore be a good approximation to ignore the error of Q in 
the analyses. 

In the following I do not distinguish between the recoil spectrum dR/ dQ and the actual 
differential counting rate dN/dQ. Since dR/dQ is usually given per unit detector mass 
and unit time, the two quantities differ only by a multiplicative constant. This constant 
can be canceled in Eq.(4.8), since it will also appear in the normalization constant H in 
Eq.(4.9). 

4.2.1 Exponential ansatz for dR/dQ 

I divide the total energy range into B bins with central points Q„ and widths 6„, 
n = 1, 2, • • • , S. In each bin, signal events will be recorded. The simulated data 
set can therefore be described by 
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The squared statistical error on dR/dQ is accordingly 
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5,000 events on Ge, 7 bins 




Q [keV] 



Figure 4.1: The curve shows the theoretical predicted recoil energy spectrum for the 
shifted MaxweUian velocity distribution fi^shi^jVe) in Eq.(3.29) with the Woods-Saxon 
form factor F^^{Q) in Eq.(3.17). The data points with error bars show simulated ex- 
perimental data produced from this spectrum (5,000 total events, = 100 GeV/c^, 
= 70.6 GeV/c^ for ^^Ge, Vq = 220 km/s, Ve = 231 km/s, and the Galactic escape 
velocity Vesc = 700 km/s as the cut-off velocity of the velocity distribution in Eq.(3.12)). 
The vertical error bars show the statistical uncertainties of the measurements, while the 
horizontal error bars indicate the bin widths. 

As mentioned at the end of the previous section, the predicted recoil spectrum re- 
sembles a falling exponential (see Eqs.(3.21') and (3.31') in Sec. 4.1). This is confirmed 
in Fig. 4.1, which shows the predicted recoil spectrum of a 100 GeV/c^ WIMP on ''^Ge 
by means of the shifted Maxwellian velocity distribution /i,sh(^^, "^e) in Eq.(3.29) and the 
Woods-Saxon form factor F^g{Q) in Eq.(3.17). This figure also shows the result of a 
simulated experiment, where the exposure time and cross section are chosen such that 
the expected number of events is 5,000; these have been collected in seven bins in recoil 
energy. Note that, in practice the velocity distribution in Eq.(3.12) should be cut off at a 
velocity fesc, since WIMPs with v > Vesc can escape the gravitational well of our Galaxy. 
The cut-off velocity has been taken to be fesc = 700 km/s. 

While an approximately exponential function can be approximated by a linear ansatz 
only over a narrow range, i.e., for small bin widths, the logarithm of this function can be 
approximated by a linear ansatz for much wider bins (some detailed discussions about 
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linear approximations can be found in App. D.4). This corresponds to the ansatz 

Here r„ is the standard estimator for dR/dQ ai Q = Qn given in Eq.(4.15), rVi is the real 
value of the recoil spectrum at the point Q = Qn, 

^"^(^)..' 

and kn is the logarithmic slope of the recoil spectrum in the n-th Q-bin. 

Now our task is to find estimators for f„ and A;„ in Eq.(4.18) using (simulated) data. 
Note that, for kn 7^ 0, f„ 7^ r„ = Nn/bn and cannot be estimated from the number of 
events Nn in the n-th bin alone. Instead, from the first part of Eq.(4.18), one has 
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where, for simplicity, I have introduced the dimensionless quantities 
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Hence, it can be found that 
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depends on kn- Moreover, using the first and second moments of the recoil spectrum in 

the n-th bin, one can find that 
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where ^^\n denotes the average value in the n-th bin. or, equivalently, kn can not 
be solved analytically by only using Eq.(4.23). They can, however, be found numerically 
once 

Q - Qnln = ^ E iQn,i ' Qn) (4-25) 

is known. Alternatively, multiplying both sides of Eq.(4.23) with bn/nn and adding to 
Eq.(4.23), one can calculate the logarithmic slopes as 



kn = , (4.26) 
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where 



{Q Qn)'^\n {Qn,i Q 



(4.27) 



n i=i 



is now estimated from the data directly. Note that kn determined either from Eq.(4.23) 
or from Eq.(4.26) is independent of the normahzation r„ or r„. 

On the other hand, the second, equivalent expression in Eq.(4.18) means that we can 
still use the quantities r„ = N^/hn as normalization. However, the logarithmic slopes, 
kn, solved by either Eq.(4.23) or Eq.(4.26) describes the spectrum dR/dQ at the shifted 
points Qs,n- Equivalence of the two expressions in Eq.(4.18) implies 
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Note that, while Qn is simply the midpoint of the n-th Q-hva. and can thus be chosen at 
will, Qs,n here is a derived quantity and depends on the logarithmic slope kn- However, the 
second expression in Eq.(4.18) combined with „ in Eq.(4.28) has two advantages. First, 
the prefactor r„ can be computed more easily than f„ in Eq.(4.22). Second, according to 
a detailed analysis [96], it has been found that, for a given bin width, one can minimize 
the leading systematic error by interpreting the estimator of kn as logarithmic slope of 
the recoil spectrum, not at the center of the bin Qn-, but at the shifted point Qs^n- Note 
that Qs^n itself depends on kn- However, this does not introduce any additional error, if 
we simply interpret Eq.(4.28) as an - admittedly somewhat complicated - prescription for 
the determination of the Q- values where we wish to estimate the logarithmic slope of the 
recoil spectrum. 

In the rest of this section I use only Q — Qn\n from Eq.(4.23) to estimate the logarith- 
mic slope kn, since it simplifies the error analysis somewhat. Note that, for using both 
Q — Qn\n and [Q — QnY\n from Eq.(4.26) to estimate the statistical errors of them 
are correlated. ^ Using standard error propagation, we have 
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The first factor above can be calculated straightforwardly from Eq.(4.23) as 
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where I have defined the auxiliary function 
g{x) = l ( 
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The error on the average energy transfer [Q — Qnin) in Eq.(4.29) can be estimated 
directly from the data, using 

1 



(7^ (Q - Qn|n) 



(Q Qn)^|n Q Qn 



(4.32) 



^In contrast, I will use Eq.(4.26) in Subsec. 5.2.2 due to some other reasons. 
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4.2.2 Windowing the data set 

From two naive linear approximations discussed in App. D.4, one can find an important 
observation that the statistical error of both estimators of the slope of the recoil spectrum 
given in Eqs.(D.34) and (D.38) scale hke the bin width to the power —1.5 (see Eqs.(D.39) 
and (D.40)). This can intuitively be understood from the argument that the variation of 
dR/dQ will be larger for larger bins (see Fig. 4.1). Moreover, a detailed analysis for the 
relation between the statistical error of kn given in Eq.(4.29) and the bin width 6„ [96] 
shows that, for small bins, the expected error again scales hke and if the bin width 
is large, the statistical error decreases even faster with increasing bin width. One would 
therefore naively conclude that the errors of the estimated slopes would be minimized by 
choosing very large bins. 

However, as mentioned above and discussed in Ref. [96], neither a linear approximation 
of the recoil spectrum nor the linear ansatz of the logarithm of the spectrum can describe 
the real (but as yet unknown) recoil spectrum exactly. The neglected terms of higher 
powers of Q — Qn will certainly induce some uncontrolled systematic errors which increase 
quickly with increasing bin width 6„. 

Using large bins has a second, obvious disadvantage: the number of bins scales in- 
versely with their size, i.e., by using large bins we would be able to estimate fi{v) only 
at a small number of velocities. Fortunately, this can be alleviated by using overlapping 
bins, or, equivalently, by combining several relatively small bins into overlapping "win- 
dows". This means that a given data point Q„ j may well contribute to several different 
windows, and hence to the estimate of fi{v) at several values of v. This can increase 
the total amount of information about fi{v) since the only information we use about the 
data points in a given window is encoded in the average recoil energy in this window 
(through the estimating of kn by Eq.(4.23)). This averaging effectively destroys informa- 
tion. By letting a given data point contribute to several overlapping windows, this loss 
of information can be reduced. 

One other obvious disadvantage of using large bins or windows is that it would lead 
to a quite large minimal value of v where fi{v) can be reconstructed, simply because the 
central value Qi, and also the shifted point Qs,i, of a large first bin would be quite large. 
This can be again be alleviated by first collecting our data in relatively small bins, and 
then combining varying numbers of bins into overlapping windows. In particular, the 
first window would be identical with the first bin. 

A final consideration concerns the size of the bins. Choosing fixed bin sizes, and 
therefore also mostly fixed window sizes, would lead to errors on the estimated logarithmic 
slopes, and hence also on the estimates of fi (v) , which increase quickly with increasing Q 
or V. This is due to the essentially exponential form of the recoil spectrum, which would 
lead to a quickly falling number of events in equal-sized bins. A try-error analysis shows 
that the errors are roughly equal in all bins if the bin widths increase linearly (some 
different variations of binning of the data set will be given in App. D.l). 
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These considerations motivate the following set-up for the mock experimental analysis. 
One starts by binning the data, as in Eq.(4.14), where the bin widths satisfy 



bn = bi + {n- 1)6 , 



(4.33) 



I.e., 



Qn Qmin ~l~ 





'{n-lf 








2 



Here the increment 6 satisfies 
2 



B{B 



(^Quiax Qmin Bbij , 



(4.34) 



(4.35) 



B being the total number of bins, and Qmax,min being the (kinematical or instrumental) 
extrema of the recoil energy. Then 1 collect up to riw bins into a window, with smaller 
windows at the borders of the range of Q. In the rest of this section and the next chapter 
1 use Latin indices n, m, • • • to label bins, and Greek indices /x, z^, • ■ ■ to label windows. 
For 1 < /J, < nw-i the /x-th window simply consists of the first [i bins; for nw ^ 1^ ^ B, the 
//-th window consists of bins /i — nM/ + 1, /x — niy + 2, /x; and for S < < -B + n^y — 1, 
the fjL-th window consists of the last nw — n -\- B bins. This can also be described by 
introducing the indices and rifj^-i^ which label the first and last bin contributing to the 
/i-th window, with 



n,. 



1, ji < nw 

H — nw + 1, A* > nw 



(4.36a) 



and 




(4.36b) 



The total number of windows defined through Eqs. (4.36a) and (4.36b) is evidently W — 
B + nw — 1, i.e., 1 < /i < 5 + nw — 1- 

As shown in the previous subsection, the basic observables needed for the reconstruc- 
tion of fi{v) in Eq.(4.8) are the number of events in the n-th bin, A^„, as well as the 
averages Q — Qn\n in Eq.(4.25). Once Nn and Q — Qn\n can be obtained, we can then 
use Eqs. (4. 15), (4.23) and (4.28) to get r„, kn and Qs^n as well as Eq.(4.29) to get the 
statistical error of k^. One can easily calculate the number of events per window as 



(4.37) 



as well as the averages 
1 



N„ 



(4.38) 
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where is the central point of the //-th window. 

One drawback of the use of overlapping windows in the analysis is that the observables 
defined in Eqs.(4.37) and (4.38) are all correlated (for nw 7^ 1)- The slope in the /i-ih 
window, kfj,, will again be calculated as in Eq.(4.23) with "bin" quantities replaced by 
"window" quantities. We thus need the covariance matrix for the Q — Q^lfj., which follows 
directly from the definition in Eq.(4.38) (detailed deviations for the covariances in this 
subsection will be given in App. E.1.1): 

cov((5-(5^|^, Q-Q^l^) 



TV TV 



r 

E (gin - gu) (gin - g|.) + Ny (g - g„i„) 



(4.39) 



where cr^ (Q — gn|n) is defined as in Eq.(4.32). In Eq.(4.39) I have assumed jJi < the 
covariance matrix is, of course, symmetric. Moreover, the sum is understood to vanish if 
the two windows //, v do not overlap, i.e., if < rii,^. 

The ansatz in Eq.(4.18) is now assumed to hold over an entire window. We again can 
estimate the prefactor as 



N,. 



w,. 



(4.40) 



W/j, being the width of the /i-th window. This implies 

1 



(4.41) 



where I have again taken ix <y. Finally, the mixed covariance matrix is given by 



cov (r^,g - g,.|i.) 



E (gu - Q\u) 



(4.42) 



Note that this sub-matrix is not symmetric under the exchange of ji and v. In the 
definition of n_ and n+ we therefore have to distinguish two cases: 

If /i < z/ : n_ = njy_, n+ = ; 

\i IJL> v : ri- = Ufj,-, n+ = riu^ . (4.43) 

As before, the sum in Eq.(4.42) is understood to vanish if n_ > n+. 

The covariance matrices involving the estimators of the logarithmic slopes A;^, derived 
from Eq.(4.23) with n ^ fi everywhere, can be calculated in terms of the covariance 
matrices in Eqs.(4.39) and (4.42): 



and 



cov {kf^, K) 



cov (r^, K) 



cov 



(g Q^i\ij,-i g gi/|i/^ ) 



cov {r^,Q- g^l^) , 



(4.44) 



(4.45) 



where is as in Eq.(4.21) with n — >■ /i, and the function g[x) has been defined in 
Eq.(4.31). 
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4.2.3 Reconstructing the velocity distribution 



Now we are ready to put all pieces together to reconstruct the velocity distribution and 
compute its statistical error. Substituting the ansatz in Eq.(4.18) (with the replacement 
n — > //) into Eq.(4.8), the reconstructed normalized one-dimensional velocity distribution 
function can be expressed as 



for /J, — 1, 2, ■ ■ 



2Q 



B + nw — 



_d_ 

dQ 



In F'iQ) 



(4.46) 



1. Here Qs^fj, is given by Eq.(4.28) with n — > /i, and, 



a 



(4.47) 



as in Eq.(3.10). Finally, the normalization constant M defined in Eq.(4.9) can be esti- 
mated directly from the data: 

-1 



2 
a 



E 



1 



(4.48) 



where the sum runs over all events in the sample. 

Since neighboring windows overlap, the estimates of fi{v) at adjacent values of are 
correlated. This is described by the covariance matrix 

fl,r{Vs,f,)flAVs. 



cov 



Til, r„ 



+ {2Ny 



cov {kfj,, ku) 



F'{Qs,u) 



cov (r^, k^) + — ^ z/) I , (4.49) 



where the covariance matrices involving the normalized counting rates and logarith- 
mic slopes k^ have been given in Eqs.(4.41), (4.44), and (4.45). In principle Eq.(4.49) 
should also include contributions involving the statistical error of the estimator for J\f in 
Eq.(4.48). However, this error and its correlations with the errors of the and k^ has 
been found to be neghgible compared to the errors included in Eq.(4.49). 

Figs. 4.2 are results for the reconstructed velocity distribution, for "typical" simulated 
experiments with 500 (top) and 5,000 (bottom) events. In the top frame B — 5 bins has 
been chosen, the first bin having a width bi = 8 kcV, and up to three bins have been 
combined into a window. Since the last bin is in fact empty, this leaves us with W = 6 
windows, i.e., we can determine /i for six discrete values of the WIMP velocity v; recall 
that these "measurements" of fi^r are correlated, as indicated by the horizontal bars in 
the figure. In the lower frame B = 10 bins with bi = 10 keV have been chosen, and up to 
four bins have been combined into one window. The bins are thus significantly smaller 
than in the upper frame. As a result, the last two bins are now (almost) empty, leaving 
us with W = 11 windows. Figs. 4.2 indicate that one will need at least a few hundred 
events for a meaningful direct reconstruction of fi{v). 
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500 events, 5 bins, up to 3 bins per window 
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5,000 events, 10 bins, up to 4 bins per window 
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Figure 4.2: The WIMP velocity distribution reconstructed from a "typical" experiment 
with 500 (top) and 5,000 (bottom) events. The smooth curves show the input distribu- 
tions, which are based on Eq.(3.29). The vertical error bars show the square roots of the 
diagonal entries of the covariance matrix given in Eq.(4.49); the horizontal bars show the 
size of the window used in deriving the given value of /i The overlap of these horizontal 
bars thus shows the range over which the values of fi^r are correlated. Parameters as in 
Fig. 4.1. 
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Furthermore, a x} distributions has been defined by 

Xf = 777 II^M^ IflAVs,!^) - flMVs,^^)] [flAVs,^) - fl,th{Vs,u)] ■ (4.50) 

Here fi^r is the estimate in Eq.(4.46) of the velocity distribution, /i^th is a theoretical 
predicted velocity distribution (e.g., the input distributions in Figs. 4.2), and C is the 
inverse of the covariance matrix of Eq.(4.49). This xj distribution allows statistically 
meaningful tests of the predicted velocity distribution function. 

More details about this xj distribution and some applications can be found in Ref. [96]. 

4.2.4 Determining moments of the velocity distribution 

As mentioned in the previous subsection, a direct reconstruction of the WIMP ve- 
locity distribution fi{v) will only be possible once several hundred nuclear recoil events 
have been collected. This is a tall order, given that not a single such event has so far 
been detected (barring the possible DAMA observation and a few candidate signals, see 
Sees. 3.7 to 3.9). The basic reason for the large required event sample is that, fi{v) 
being a normalized distribution, only information on the shape of fi{v) is meaningful. In 
order to obtain such shape information via direct reconstruction, we have to separate the 
events into several bins or windows. Moreover, each window should contain sufficiently 
many events to allow an estimate of the slope of the recoil spectrum in this window. 

On the other hand, at the end of Sec. 4.1 1 have also given expressions for the moments 
of fi{v) in Eqs.(4.11) to (4.13). With the exception of the moment with n — —1, these 
are entirely inclusive quantities, i.e., each moment is sensitive to the entire data set; no 
binning is required, nor do we need to determine any slope (with one possible minor 
exception; see below). It thus seems reasonable to expect that one can obtain meaningful 
information about these moments with fewer events. 

The experimental implementation of Eq.(4.11) is quite straightforward. For Qthre — 0, 
the normalization Af has already been given in Eq.(4.48). The case of non- vanishing 
threshold energy Qthrc can be treated straightforwardly, using Eq.(4.13). To that end 
one needs to estimate the recoil spectrum at the threshold energy. One possibility would 
be to choose an artificially high value of Qthre? and simply count the events in a bin 
centered on Qthrc. However, in this case the events with Q < Qthre would be left out of the 
determination of the moments. We therefore should keep Qthrc as small as experimentally 
possible, and to estimate the counting rate at threshold using the ansatz in Eq.(4.18). 
Since we need the recoil spectrum only at this single point, we only have to determine the 
quantities ri and ki parameterizing dR/dQ in the first bin; this can be done as described 
in Subsec. 4.2.1, without the need to distinguish between bins and "windows". Introduce 
the shorthand notation 

rthre ^ (^] = rie^^W— «-^) . (4.51) 
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Then, combining Eqs.(4.11) and (4.13), the n-th moment of the velocity distribution 
function can be rewritten as 

-1 



thrc J 



9^(n+l)/2 
^Vthrc 



thre J 



+ {n + 1)4 



(4.52) 



where the integral /„ defined in Eq.(4.12) can be estimated through the sum: 

r)(n-l)/2 

as Eq.(4.48). Since all /„ are determined from the same data, they are correlated with 

g(n+m-2)/2 



COv(/„, Im) = Yl 



F\Qa 



(4.54) 



This can e.g., be seen by writing Eq.(4.53) as a sum over narrow bins, such that the 
recoil spectrum within each bin can be approximated by a constant. Each term in the 
sum would then have to be multiplied with the number of events in this bin; Eq.(4.54) then 
follows from standard error propagation. Note that, when re-converted into an integral, 
the expression for cov(/o, /q) will diverge logarithmically for Qthre 0. Equivalently, 
the numerical estimate of this entry can become very large if the sample contains events 
with very small Q- values. But, according to some numerical simulations, there should be 
no problem for samples with Qthre > 1 keV. Many existing experiments in fact require 
significantly larger energy transfers in their definition of a WIMP signal. 

In order to calculate the statistical error of {v'^) in Eq.(4.52), one needs at first the 
error of rthre which can be obtained from Eq.(4.51) as 



o-^(rthr 



' thre 



+ 



Qthre ~ Qs,l — ki 



'dQs 



n 2 



dki 



(4.55) 



Here the squared errors for ri and ki are simply the corresponding diagonal entries of the 
respective covariance matrices given in Eqs.(4.41) and (4.44), and the definition of Qg^i 
in Eq.(4.28) implies 

^dQsA ^ 1 ^ fbi\ fhki^ 



QsA + ki 



dki 



(4.56) 



where Qi is the central Q-value in the first bin. It should be noted that the first term 
in Eq.(4.11) is negligible for all n > 1 if Qthrc — 1 keV. However, even for this low 
threshold energy it contributes significantly to the normalization constant J\f, as described 
by Eq.(4.13). Of course, the first term in Eq.(4.11) always dominates for n = —1. This 
is not surprising, since the very starting point of this analysis, Eq.(3.12), already shows 
that the counting rate at Qthre is proportional to the "minus first" moment of the velocity 
distribution. 

One needs also the correlation between the errors on the estimate of the recoil spec- 
trum ai Q = Qthre and the integrals /„. It is clear that these quantities are correlated. 
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since the former is estimated from all events in the first bin, which of course also con- 
tribute to the latter. These correlations can be estimated by using the ansatz in Eq.(4.18), 
which makes the following prediction for the contribution of the first bin to the integrals: 



^,l = n 



This immediately implies 

dln,i _ In,l 
dri ri ' 

and 



^kiiQ-Qs,i) 



dQ. 



dln,l 



Qs,i + ki 



'dQ 



(4.57) 



(4.58a) 



(4.58b) 



Note that /„ i and In+2,1 in Eqs. (4.58a) and (4.58b) can be evaluated as in Eq.(4.53), 
with the sum extending only over events in the first bin: 



4,1 - 



The correlation between rthre and /„ is then given by 

COv(rthrc, In) 



(4.59) 



^thre In,l 



rl 



+ 



Qthre ~ Qs,\ ~ ki 



'dQ 



X 



'n+2,1 

In,l 



Qs.l ~ h 



'dQs 



dh 



a^h) . 



(4.60) 



Finally, these ingredients allow us to compute the covariance matrix for the estimates 
of the moments of the velocity distribution: 

cov((i;"),(i;™)) 
= [(^^") (^"^>cov(7o, lo) + a"+"^(n + l)(m + l)cov(7„, Im) 

- a^{m + l)(t;")cov(/o, /„) - a'\n + l)(t;™)cov(/o, 4) 

+ D^D^a^ (rthre) - (i^m(^") + £'n(^'"))cOv(rthre, h) 

+ Q;'"(m + l)D„COv(rthre, Im) + Q;"(n + 1)D^C0V (rthre, In) ■ (4.61) 

Here I have introduced the modified normalization constant: 

A4r^(f)A/', (4.62) 

which exploits the partial cancellation of the a factors between Eqs. (4. 11) and (4.13), 
and the quantities 



'div^^y 



AAm \drthrej -^^(Qthre) 



a 
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(4.63) 



Note that, in practice, one can determine {v'^) by a single experiment with a large number 
of events, or by averaging over many experiments with a relatively small number of events. 
However, numerical simulations [96] show that in the second case the average values of the 
reconstructed moments do not exactly converge to the input (exact) values. In order to 
understand this, consider the simple case Qthrc = 0. The moments arc then proportional 
to the ratio In/Io (see Eq.(4.52)). The distortion arises because {In/Io) {^n)/{Io), 
where the averaging is over many simulated experiments. The leading correction terms 
for small Qthrc and not very large first bin can be found as (a detailed calculation by 
using Taylor expansion to second order will be given in App. E.2) 



= a'^AA2|(n + l)[cov(/o,/„) - A4./„cov(/o, Jo) 



+ 2 



^(n+l)/2 
Vthrc 



COv(rthre, Iq) " nhreA/mCOv(/o, /q) 



(4.64) 



where the second line in Eq.(4.60) is significant only for n — —1. Note that this correction 
becomes very small if the statistical errors on the /„ as well as on rthre become small. 

Meanwhile, according to some detailed numerical analyses [96], an "error on the error" 
should be added. The contribution to the diagonal entries of the covariance matrix given 
in Eq.(4.55) can be estimated as 



(cov(/„, /„)) = Yl ■ 



(4.65) 



the off-diagonal entries are then scaled up such that the correlation matrix remains unal- 
tered. The numerical analyses show also that very rare events with large recoil energies 
contribute significantly more to the higher moments. Hence, an experiment with a small 
number of events will usually underestimate {v'") and, especially, its error; the problem 
will become worse with increasing n. However, because this method uses whole exper- 
imental data together to determine the moments of fi{v), it has also been found that, 
based only on the first two or three moments, some non-trivial information can already 
be extracted from 0(20) events. 

More details and discussions about the reconstruction of the velocity distribution and 
determination of its moments can be found in Ref . [96] . 



4.3 Determining the WIMP mass 

In the previous two sections I discussed how to use a recoil spectrum from direct Dark 
Matter detection as well as experimental data directly (i.e., the measured recoil energies) 
to reconstruct the velocity distribution function of WIMPs as well as to determine its 
moments. As noted earlier, for both of these reconstruction methods we need to know 
the mass of the incident WIMPs m^. In well-motivated WIMP models from elementary 
particle physics, can be determined with high accuracy from future collider experiment 
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data. However, one has to check experimentally that the particles produced at colliders 
are in fact the same ones seen in Dark Matter detection experiments which form the 
Galactic halo. In this section I present a method for (self-)determining the WIMP mass 
based on the determination of the moments of the velocity distribution function, (v"), 
(presented in Sec. 4.2.4) from two (or more) experimental data sets with different target 
materials. ^ 



4.3.1 Neglecting Qthve 

As mentioned in the end of Sec. 4.1, the basic idea for using two different detector 
materials to determine the WIMP mass is that, from independent direct WIMP detection 
experiments with different target nuclei, the measured recoil spectra should lead to the 
same (moments of the) velocity distribution function of incident WIMPs. 

For the case that the threshold energy Qthre can be neglected, the n-th moment of the 
velocity distribution function, {v"'), in Eq.(4.52) can be expressed simply as 

{v^) = a'\n + l) (y^) , (4.66) 

where /„ and Iq can be estimated by Eq.(4.53). Suppose X and Y are two target nuclei. 
Eq.(4.66) implies 

Note that the form factor F'^{Q) in Eq.(4.53) for estimating /„x and /„ y are different. 
Then, according to the definition of a in Eq.(3.11) with the expression of the reduced 
mass mr in Eq.(3.6) and using some simple algebra, one can find the WIMP mass as 

= , , (4.68) 

"^n - \Jmx/mY 

where I have defined 

7^„^^=f^•f^V'^ n^O, -1. (4.69) 

Fig. 4.3 shows the ratios of the reproduced WIMP masses estimated by Eq.(4.68) with 
different combinations of target nuclei to the input (true) one as functions of the input 
WIMP mass. ■^^Si, ^°Ar, and ^^Ge have been chosen as three target nuclei and thus three 
combinations for defined in Eq.(4.69) with n — 1 are shown. TZn has been estimated 
by the integral form of /„ in Eq.(4.12) with a maximal measuring energy of 200 keV. 
The theoretical predicted recoil spectrum for the shifted Maxwellian velocity distribution 



"^Note that the ansatz here is quite different from that used in Ref. [97], which assumes different 
WIMP velocity distributions with two input parameters: the WIMP mass and the WIMP-nucleon cross 
section, and then analyses with which precision in the usual WIMP mass-cross section plane the WIMP 
mass can be reproduced from the direct detection experiment. 
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Figure 4.3: The curves show the ratios of the reproduced WIMP masses estimated by 
Eq.(4.68) with different combinations of target nuclei to the input (true) one as functions 
of the input WIMP mass. with n = 1 has been estimated by the integral form 
of In with a maximal measuring energy of 200 keV. The recoil energy spectrum for a 
shifted Maxwellian velocity distribution with the Woods-Saxon form factor has been 
used (parameters as in Fig. 4.1). The solid (blue) line, the dashed (black) line, and the 
dash-dotted (red) line are for ^°Ar + ^^Si, ^^Ge + ^°Ar, and ™Ge + ^^Si combination, 
respectively. The straight dash-dotted (green) line denotes 1. 

function, {dR/dQ)sh in Eq.(3.31), with the Woods-Saxon form factor -Fws(Q) Eq.(3.17) 
has been used. In Fig. 4.3 one can see obviously a deviation of the reproduced WIMP mass 
from the input (true) one as input > 60 GeV/c^. The heavier the nuclear masses 
of two target nuclei, e.g., ^^Ge -|- ^'^Si, the larger the deviation from the true WIMP 
mass. This is caused by introducing the maximal measuring energy for estimating J„. As 
discussed in Subsec. 4.2.4, the heavier the nuclear mass mN, or, equivalently, the larger 
a, and the larger n, the more the contribution to J„ comes from the high Q region, and, 
for a fixed maximal measuring energy, the smaller the value for and then for will 
be estimated. As shown in Fig. 4.3, for n = 1 and input = 200 GeV/c^, the deviation 
with Qmax = 200 keV is around 20%. However, according to the numerical analysis, 
with Qmax = 250 keV or 300 keV, this deviation will be reduced to around 10% or even 
only 5%. Later we will see, due to a quite large statistical error with very few events, a 
deviation around 10% for input = 200 GeV/c^ is not very bad. Moreover, for input 
< 120 GeV/c^, the deviation should be less than 5% or even 1% for Qmax = 200 keV 
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□max = 200 keV, Qmin = 1 keV, n = 1 , 25 + 25 events 




Figure 4.4: The curves show the statistical errors estimated by Eq.(4.70) with different 
combinations of target nuclei as functions of the input WIMP mass. Each experiment has 
25 events, i.e., totally 50 events. Parameters and indications of the lines as in Fig. 4.3. 



or 250 keV. 

Furthermore, the statistical error on the reproduced WIMP mass can be obtained 
from Eq.(4.68) directly as 
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1/2 



(4.70) 



where <J^ {In,x) = cov (/„ -^n,x) and cov {Iq x, In,x) and so on can be estimated from 
Eq.(4.54). 

Fig. 4.4 shows the statistical errors estimated by Eq.(4.70) with three different com- 
binations of target nuclei as functions of the input (true) WIMP mass. Each experiment 
has 25 events, i.e., totally 50 events. Note that, in order to use the integral form of 
cov(/„, Im) in Eq.(4.54), a threshold energy Qmin = 1 keV has been given. In Fig. 4.4 one 
can observe that the larger the mass difference between two detector nuclei, the smaller 
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the statistical error will be. Hence, the combination with the largest mass difference, e.g., 
^^Ge + ^^Si will have the smallest statistical error. In principle one other combination: 
^^^Xe + ^'^Ar has larger mass difference and should have an even smaller statistical error. 
However, because the Woods-Saxon form factor has been used here, the integral form of 
cov(/„, Im) in Eq.(4.54) has a pole at Q ~ 100 keV! Thus ^^^Xe has been not used for this 
simulation. On the other hand, despite of the factor l/\n\ in Eq.(4.70), it has been found 
that the statistical errors increase with increasing n, except the ^°Ar + ^*Si combination. 
For this combination, the statistical error with n = 2 is a little smaller than with n = 1; 
but, with n = 3, the statistical error is significantly larger (and, as discussed above, the 
deviation of the reproduced WIMP mass should also be larger). Hence, n = 1 should be 
the best choice for and crlm^) in Eqs.(4.68) and (4.70), respectively. 

Figs. 4.5 show the reproduced WIMP mass with the statistical error by using ''^Ge 
and ^^Si as two target nuclei as a function of the input (true) WIMP mass. From the 
upper frame, it can be found that, despite of the very few (25 + 25, totally 50) events and 
correspondingly very large statistical error, for < 100 GeV/c^, one can already extract 
some meaningful information on the WIMP mass. For example, for — 25 GeV/c^ and 

= 50 GeV/c^ we will reproduce ~ (25± 13) GeV/c^ and ~ (50 ±31) GeV/c^. 
For the case with 500 (250 + 250) total events, the statistical error will be reduced to less 
than 5 and 10 GeV/c^, respectively! Certainly, as shown in the lower frame of Figs. 4.5, 
for the case with 500 total events, the deviation of the reproduced WIMP mass from the 
input one becomes important. Nevertheless, in practice, an experiment with more than 
200 events should have a larger maximal measuring energy, and, as discussed above, the 
deviation can (should) be strongly reduced. 

For the simplified simulations with the integral form of /„ presented above, the event 
numbers from both experiments have been considered to be equal. Practically, as de- 
scribed in Subsecs. 3.4.1 to 3.4.3, experiments with the higher mass nuclei, e.g., Ge or 
Xe, are expected to measure (much) more signal events. However, according to the ex- 
pression for o"(m-^) in Eq.(4.70) and the definition of TZn in Eq.(4.69), it can be found 
that only the terms in the brackets depends on the event number and the contributions 
from the two experiments are independent of each other. Moreover, a detailed analysis 
of contributions from different terms of cr(m;^) shows that the prefactor 



which depends practically only on the choice of the combination of the two target nuclei 
is very large for every combination, while the terms in the brackets with the factor l/\n\ 
are actually quite small. This implies that one can not reduce the statistical error of 
estimated by Eq.(4.70) by improving only one experiment with even very large event 
number, since the contribution from the other (poor) experiment will dominate the error. 




(4.71) 
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Qmax = 200 keV, Qmin = 1 keV, n = 1 , 25 + 25 events, Ge-76 + Si-28 




Qma>! = 200 keV, Qmin = 1 keV, n = 1 , 250 + 250 events, Ge-76 + Si-28 




50 100 150 

inchi_in[GeV] 



Figure 4.5: The reproduced WIMP mass with the statistical error by using Ge and Si 
as two target nuclei as a function of the input WIMP mass. The solid (red) line indicates 
the reproduced WIMP mass estimated by Eq.(4.68), the dashed (red) lines indicate the 1- 
0" statistical error estimated by Eq.(4.70). The straight dash-dotted (green) line indicates 
the input (true) WIMP mass. Each experiment has 25 (250) events, i.e., totally 50 (500) 
events, in the upper (lower) frame. Parameters as in Fig. 4.3. 
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4.3.2 With Qthre > 



For the case that Qthre in Eq.(4.52) can not be neglected, TZn defined in Eq.(4.69) 
should be modified to the following general form: 

"2gtt!i^Vthre,X + {n+ l)/n,xFl(Qthre,x)' 



T^niQ thre) 



2<5jhre,X^thre,X + -^0,X-^l(Qthre,x) 

2Qthre thre,y + -^0,F-^y (Qthre.y) 



X 



2Qihrey ^thre,y + (n + l) -^n,y-Fy (Qthre.y) 



l/n 



(4.72) 



where rthrc,x and rthrc.y should be determined by Eq.(4.51) (practically) with different 
ri, ki, Qs,i, and Qthre- In this general form of 7^n(Qthre) there are totally six variables: 
i^thrc,x, In,x, -^o,x and the other three for nucleus Y. This should generally produce a 
larger statistical error than that estimated by Eq.(4.70) due to the contribution from 
rthrc (the statistical error of TZniQthrc) will be given in App. E.3). However, one can 
practically reduce the number of variables by choosing n = —1: 

20thre,X^thre,X + -^0,X-^l(Qthre,x) 



'^-l(Qthre) 



r thrc.y 



^'thre.X _ 2(5^/^g yrthre,y + -^0,y-^y (Qthre.Y ) _ 

Then a{TZn) in the first fine of Eq.(4.70) should be replaced by 

(7(7^_l((5thre)) 



(4.73) 



= 7^-l((5thre)' 



1 2 



-2<5tYre,x'^thre,X + -^0,X-^l(<5thre,x) - 



X 



t^Mnhre,x) , Cr'^{h,x) 2C0V (rthre,X, -^0,x) 



J.2 

' thre,X 



+ 



-'o,x 

^ 1/2 



'"thre,X-^0,X 



(4.74) 



where o"^ (rthre.x) and cov (rthre,x, -^o,x) and so on can be calculated by Eqs.(4.55) and 
(4.60). 
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Chapter 5 

Annual Modulated Event Rate 



In the previous chapter I have presented methods for reconstructing the velocity 
distribution function and its moments from the time-averaged recoil energy spectrum 
fitted to experimental data as well as from data directly. The annual modulation of the 
event rate discussed in Sec. 3.2 has been ignored. As shown in Sec. 4.2, in the foreseeable 
future with rare signal events, the statistical errors will remain large and thus this is a 
reasonable first approximation. However, for the future detectors with strongly improved 
sensitivity and (very) large target mass (large exposure) , the formulae and methods have 
to be extended to allow for an annual modulation of the event rate. 

In the first section of this chapter I extend the method developed in the previous 
chapter by considering an arbitrary, but cosine-like time- dependent recoil spectrum with 
a one-year period. In the second section I present the method for reconstructing the 
amplitude of the (possible) annual modulation of the velocity distribution function. An 
alternative, better way for checking the annual modulation of the event rate will also be 
described. 



5.1 Taking into account the annual modulation 

For simplicity, in this chapter I take = in Eq.(3.30) and rewrite it as 

v^{t) = i;o[l.05 + 0.07cos(a;i)] , (3.30') 

with 

This means that in the following analyses experiments (data) have been assumed to 
be started (collected), i.e., t = 0, when Ve is maximal (around June 2nd, theoretically 
predicted) and the time t will be measured in unit of "day" . 

As discussed in Sees. 3.1 and 3.2, roughly speaking, the differential event rate for 
direct WIMP detection is proportional to the WIMP flux, or, equivalently, the velocity of 
the Earth relative to the WIMP halo. And, due to the motion of the Earth on an elliptical 
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orbit around the Sun, the projection of the Earth's orbital speed on the orbital speed 
of the Sun around the Galactic center is approximately a cosine function. Therefore, as 
shown in Eqs.(3.30) and (3.30') above, the differential event rate should theoretically be a 
cosinusoidal function (c.f. the DAMA results in Figs. 3.3). On the other hand, substituting 
Ve{t) in Eq.(3.30') into Eq.(3.31), it can be found that {dR/dQ)sh are not exact cosine 
but cosine-like functions with a period of 365 days (shown in Figs. 5.1). According to 
this observation, I assume generally an arbitrary, but cosine-like time- dependent recoil 
spectrum with a one-year period and then expand this spectrum and its corresponding 
velocity distribution function as Fourier cosine series as: 

fdR\ fdR\ (dR\ , , (dR\ ,^ , ,^ 

and 

/i {V:t) ^ /i,(o) (v) + /i,(i) (v) cos{u;t) + /i,(2) (v) cos(2u;i) H . (5.3) 

According to Eq.(3.12), {dR/dQ)t and fi{v,t) must satisfy the equation for the time- 
dependent WIMP-nucleus scattering spectrum: 

and, consequently, each pair of their Fourier coefficients must satisfy a time-independent 
equation: 

Moreover, if we neglect (due to the very low detection rate discussed in Sec. 3.1 and the 
tiny difference shown in Figs. 5.1, we can practically neglect) the terms with m > 2 in 
Eqs.(5.2) and (5.3), {dR/dQ){Q) and fi,(p){v) above are the time-averaged scattering spec- 
trum and the time-averaged velocity distribution function of WIMPs, which we considered 
in Chap. 4, and {dR/dQ)(^i) and fi,(i){v) are the amplitudes of the annual modulations 
of the scattering spectrum and its corresponding velocity distribution. In addition, since 
{dR/dQ){rn) are Fourier coefficients of {dR/dQ)t, we have 



and 




Now, as mentioned in Subsec. 4.2.2, the important elements needed for the recon- 
struction of fi^r in Eq.(4.46) are the number of events Nj^ in the /i-th window given in 



fi{v,t) 



dv , 



(5.4) 



/i,m(^^) 



dv , 



m 



0, 1, 2, 



(5.5) 
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Figure 5.1: The solid (red) curves are the predicted modulations of the recoil energy 
spectrum for the shifted Maxwellian WIMP velocity distribution, {dR/dQ)sh in Eq.(3.31), 
with the Woods-Saxon form factor F^^{Q) in Eq.(3.21); the dashed (blue) curves are 
cosine functions with an amplitude [{dR/dQ)sh{t = 0) — {dR/dQ)sh{{t = 7r/2)]/2. Here 
I have used = 100 GeV/c^, = 70.6 GeV/c^ for ^^Ge, vq = 220 km/s. The upper 
and lower frames are drawn for Q = 15 keV and Q = 30 keV, respectively. 
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Eq.(4.37), as well as the averages Q — Q^j^ given in Eq.(4.38), which are theoretically 
defined as, respectively, 



and 

where I have used generally the A-th moment of the averaged recoil spectrum {dR/ dQ){Q). 
Substituting Eq.(5.6) into Eqs.(5.8) and (5.9), it can be found easily that, for a time- 
dependent recoil spectrum with a one-year period, 

365 Jo JQa-wa/2 \dQ L 365 



and 

{Q - Q,Y 



" N 



1 /-aes i-Q^+w^/2 
565 Jo JQu-wu 2 




N 1 



365 Jo Jq^-w^/2 



where Q/i,*, ^ = 1, 2, A^^^^i yr, are the measured recoil energies from the direct 

WIMP detection experiment in the /i-th window in one year. Note that the sign 
in the second line of Eq.(5.11) denotes not mathematically equal but an experimental 
estimator for (Q — <5/i)'^U■ 
Comparing these results with the expressions in Eqs.(4.25) and (4.27), Eqs.(5.10) and 
(5.11) show that, for an arbitrary time- dependent recoil spectrum with a one-year period 
(even though it is not cosine- like) , we just have to take the experimental data over some 
whole years to find out the average event number (per day) and the annual average value 
of the energy transfer (Q — Q^)'^ in the /i-th window (or bin). Then we can reconstruct 
the time-averaged velocity distribution by means of the method presented in the previous 
chapter directly. Moreover, the results above show also that it is actually not necessary 
to know when Ve is maximal but only use all events collected in these (whole) years. 



5.2 Reconstructing the modulated amplitude of ) 

In this section 1 follow the trick used with {dR/dQ){Q) in the previous section and 
develop a method for reconstructing the (annual) modulated amplitude of fi{v). Mean- 
while, I will also introduce two criteria for checking the annual modulation of the event 
rate. 
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5.2.1 Criteria for the annual modulation 



Replacing {dR/ dQ)(Q\ in Eq.(5.9) by {dR/dQ)n\ in Eq.(5.7), it can be found that ^ 



Qn+bn/2 



Nn JQ„-bn/2 



{Q-Qn 



'dR\ 



dQ 



(1) 



365 rQn+bn/2 



365 JO JQn-b„/2 



fdR? 



dQ dt 



1 yr 



iQn,i - Qn) COs{ujtn,i) , 



(5.12) 



^^",1 yr i=l 

where tn,i, i — I, 2, iV„,i yr, n = 1, 2, B, are the "measuring times" at 

which we measure the recoil energies Qn,i in the n-th bin in one year. Note that, first, 
the "=" sign in the second line of Eq.(5.12) denotes again an experimental estimator; 
second, the factor cos(a;t„^j) comes from the integral in Eq.(5.7) (not from the second 
term of Ve{t) in Eq.(3.30')!). On the other hand, by substituting Ve{t) in Eq.(3.30') 
into {dR/dQ)sh in Eq.(3.31), it can be found that the ratio of the modulated amplitude 
of the recoil spectrum, {dR/dQ)(^i), to the time-averaged recoil spectrum, {dR/dQ)(^Q^, 
increases monotonically with the recoil energy Q and is approximately a linear function 
of Q (shown in Fig. 5.2). Hence, I introduce an ansatz for the modulated amplitude of 
the recoil spectrum in the n-th bin: 

fdR\ fdR\ 



(0) 



^n(Q Qn) ~l~ 



n 



1, 2, 



B. (5.13) 



Note that {dR/dQ)(Q) here indicates generally a timc-avcragcd recoil spectrum, not spec- 
ified to the exponential ansatz in Eq.(4.18). Substituting the ansatz for {dR/ dQ){i)^n into 
the left-hand side of the Eq.(5.12), it can be found that 



Nn JQn-b„/2 



{Q - QnY 



'dR\ 



dQ ^IniQ- Qn)^+\ + hn{Q- Qn)''\n ■ (5.14) 



AQJil),n 

Setting A = and 1 and combining Eqs.(5.12) and (5.14), 1^ and /i„ in Eq.(5.13) can be 
solved as 



2{Q - Qn) cos(cjt) |„ - 2 cos(cjt) |„ Q - Qn 



{Q Qn)'^\n Q Qnln 



and 



hn = 



2cOs{u;t)\n{Q-Qn)^ 


n-2{Q - Q„) COs{ujt)\n Q - Qn n 


{Q- 


Qn)' 


n Q Qn 


2 
n 



where I have defined 



N, 



n,l yr 



{Q -Qn)^COSP{ujt)\n 



Nn,l 



(Qn,i - Qn) COSP{utn,i) ■ 



(5.15) 



(5.16) 



(5.17) 



yr i=l 



^Foi simplicity and clarity, I discuss in this section the case with bins, thus /x and in Eq.(5.9) have 
been replaced here by n and 6„, respectively. However, all formulae given in this section can be used for 
the case with windows by substituting n and 6„ by /z and w^. 
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0.251 



250 



Q[keV] 



Figure 5.2: The curve shows the ratio of the modulated amphtude of the recoil spectrum, 
{dR/dQ)(^i), to the time-averaged recoil spectrum, {dR/dQ)[o), as a function of the recoil 
energy Q. Parameters as in Figs. 5.1. 

Due to the annual modulation effect, around t — we should get more events than 
around t — tt. Recall that I have assumed that experiments start when is maximal. 
Thus, even though J^'^ cos{u;t) dt — 0, cos(a;i)|„ above are generally not equal to 0. 
Moreover, Z„ in Eq.(5.15) can be rewritten as 



where cov{Q — cos(a;i)|„) is the covariance between the average value of the mea- 
sured recoil energies Qn,i — Qn and that of cos(a;f„_i) in the n-th bin. According to Fig. 5.2 
and the ansatz in Eq.(5.14), Z„ should be generally > and and this means that Q — Qn\n 
and cos{ijjt)\n should be positively correlated. This result offers a better way to test the 
annual modulation effect! Traditionally, in order to confirm the annual modulation of 
the event rate, one has to collect the recoil events in a given energy range in several short 
time intervals (a few days to a couple weeks) and than compare the numbers of collected 
events in the different time intervals in one year (e.g., the DAMA 4-year and 7-year results 
shown in Figs. 3.3). As mentioned in Sec. 3.2, the annual modulation of event rate is 
expected to be only a few percent (about —4% ~ 5% in the energy range between and 
50 keV, see Fig. 5.2) and this method can be used once more than one hundred events 
(in a few days!) have been accumulated. However, according to the discussion above, one 
can now collect all recoil events in a relatively larger energy range (since the calculation 




(5.15') 



97 



with bins can be extended directly to that with windows) in one year (or even several 
years) and then only has to check the following quantities: 

At = cos(a;t) = ^ cos{ujta) , (5.18) 

^^tot a 



and 



Aq^t = cov [Q, cos{ut 



Q cos{ut) — Q cos{ut) 



(5.19) 



Ntot - 1 

where the averages are over all events in the energy range which one concerns and Ntot 
is the total event number in this energy range. If the annual modulation of event rate 
exists, one should than get Aj 7^ and Aq ^ > 0. Note that, for the case that some time- 
independent background events mixed into the true signals, the two quantities above will 
be underestimated through the averaging; or even worse, if most of the background events 
have been discriminated, then the contribution from the rest events can not cancel each 
other any more. However, for the case that the time-independent background events 
dominate the whole data set, one can use a quantity modified from At in Eq.(5.18): 

A[ = Y,cos{uta). (5.20) 

a 

The quantity A[ defined here is not the average but the sum of the cosine values of the 
measuring times. Hence, the contributions from background events will cancel each other 
and not change the value of A^ (too much). But in the statistical error of A^: ^ 

a\A't) = J2cos\ujta), (5.21) 

a 

there can not be any cancellation, i.e., (A'J will increase with the total event number 
Ntot- In addition, for a "time-dependent" background the quantity in Eq.(5.20) can not 
be used any more. 

Furthermore, for the check of the quantities in Eqs.(5.18) and (5.19), it is important 
to know when Ve is maximal. This offers in practice a possibility to determine (to check 
theoretically predicted) tp in Eq.(3.30). One sets the starting date of the experiment on 
January 1st, inserts a phase (f into Eqs.(5.18) and (5.19), and then finds out when the 
quantity 



At-^^ cosoj(t- if) (5.18') 

is (almost) equal to 0, which corresponds to ± 7r/2a;, and when the quantity 

AQ,t-<^ = cov (q, cosou{t - (pfj (5.19') 

has a maximal value (positive), a minimal value (negative), or is (almost) equal to 0, 
which correspond to tp, ± tt/o;, or tp±7r/2u!, respectively. Certainly, for the case that 
such annual modulation of the event rate does not exist, one will find that At and Ag^t 
are independent of (p and always (approximately) equal to 0. 



^Here I have assumed that Ntot ^ 1 and then used (cos(a;<)) = cos'^{Lot) ~ cos{ujt) /Ni 



Hot- 
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5.2.2 Reconstructing the modulated amplitude of fi{v) 

According to Eq.(4.8), each coefficient of tlie time-dependent velocity distribution 
function, fi{v,t), in Eq.(5.3), can be solved from Eq.(5.5) as 



/i,m(Q)=-AA<|-2Q- A 




m = 0, 1, 2, 



since, for all m, 



fi, 



(m) \ 



OO 



0, 



(5.22) 



(5.23) 



see Eq.(4.3). Substituting the ansatz for {dR/dQ)(^i^^n in Eq.(5.14) at first and then the 
ansatz for [dR/dQ)(^o),n in Eq.(4.18) into Eq.(5.22), the ratio of the modulated amplitude 
of the velocity distribution function to the time-averaged one at the point Q = Qn (not 
ai V = Vnl) can be obtained as (a detailed derivation will be given in App. E.4) 



Vn 



/l,(l),n(Qn) 
fl,{0),n{Qn) 



kn 



Q=Qn 



(5.24) 



Note that the first term in the brackets has been evaluated a,t Q — Qn (not at Q = Qs,n-)- 
The result in Eq.(5.24) has three advantages. First, since rjn here is the ratio of the 
modulated amplitude of the velocity distribution function to the time-average one in the 
n-th Q-hin, it is not necessary to combine the data from all bins to get the normalization 
constant, each one of these rjn is independent of the others. ^ Second, due to the same 
reason, one can estimate the rjn even if he can not obtain (enough) events in the high 
energy range (> 100 keV). It is, in contrast, necessary to collect enough events until high 
energy (^ 200 keV for a WIMP mass ~ 100 GeV/c^ and ^^Ge as target nucleus) in order 
to determine the normalization constant A/" in Eq.(4.48). Third, the rjn are independent of 
a defined in Eq.(3.11), or, equivalently, independent of the WIMP mass m^. Certainly, 
one can use the method described in Sec. 4.3 to determine the WIMP mass. But for 
the case with very rare total events or too few events in the energy range higher than 
e.g., 100 keV, the deviation of the estimation of the WIMP mass from the true one and 
the statistical error will be very large. However, due to the independence of the rjn of 
a, the reconstruction by Eq.(5.24) will not be affected by the (large) uncertainty of the 
estimation of the WIMP mass. 

The statistical errors of rin in Eq.(5.24) can be expressed as 

^"(Vn) = E f ' Ay^,n) + E I £^ I I £^ (?/^,n, y.,n) ■ (5.25) 



i/=l 



Here I have defined 



Vv^n^Q -Qn\n, (Q-Qn)^|n, COs(a;t)|„, (Q - Qn) COs(u;t) |„ . 



(5.26) 



"^For the case with windows, the r\n are no more independent of each other. But the off-diagonal 
entries of the correlation matrix of the statistical error of the r\n are practically not important. 
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and one has 

{yi^,n) = GOV y^^n) , 

with 

GOV ((Q - Qn)^COSP{cut)\n, {Q - Qn)'^ COS^M) |„) 
1 



(5.27) 



Nn-l 



{Q-QnV+''cosf+-(u;t)\, 



{Q - QnV COSP{ujt)\n{Q - QnY COS^ {l^t)\n 



(5.28) 



According to Eqs.(4.26), (5.15), (5.16), and (5.24), the derivatives of rjn with respect to 
each of the y^^n in Eq.(5.26) can be found easily as 



dQ Qyi I n ^r, 



Q Qn I n ~l~ -^n 



I K'^ 



Q Qn\ 



dQ Qn\n 

9cos(a;i)|„ 



dVn 



d{Q - Qn) cos(c^t)|„ 
where I have defined 

C^n = {Q - Qn\n) 

and 

" d 



{Q Qn)^\n ~^ Q Qn\n^n 



(^Q Qn\n 



Nn-l 



(Q Qn)^|n {Q Qn\nj 



Q=Qu 



(5.29a) 
(5.29b) 
(5.29c) 
(5.29d) 

(5.30) 

(5.31) 



''Since /„ and /i„ arc functions of both Q — Qn\n and {Q — Qn)^|n) for simplicity, I have used here 
the expression for kn in Eq.(4.26). 
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Chapter 6 

Summary and Conclusions 



In this thesis I have presented methods which allow to extract information on the 
WIMP velocity distribution as well as on the WIMP mass from the recoil energy spectrum 
dR/dQ measured in elastic WIMP- nucleus scattering experiments. In the long term the 
information on the WIMP velocity distribution can be used to test or constrain models 
of the dark halo of our Galaxy; this information would complement the information on 
the density distribution of WIMPs, which can be derived e.g., from measurements of the 
Galactic rotation curve. Meanwhile, the information on the WIMP mass can be used 
to constrain e.g., SUSY models in the elementary particle physics and compare with 
information from future collider experiments. 

In Sec. 4.1 I have derived the expression that allow to reconstruct the normahzed one- 
dimensional velocity distribution function of WIMPs, fi{v), given an expression (e.g., a 
fit to data) for the recoil spectrum. I have also derived formulae for determining the 
moments of fi{v). All these expressions are independent of the as yet unknown WIMP 
density near the Earth as well as of the WIMP-nucleus cross section. The only information 
about the nature of WIMPs which one needs is the WIMP mass. 

Then, in Sec. 4.2, I have presented methods that allow to apply the expressions derived 
in Sec. 4.1 directly to experimental data, without the need to fit the recoil spectrum to 
a functional form. A good variable that allows direct reconstruction of fi{v) is the 
average recoil energy in a given bin (or "window", introduced in Subsec. 4.2.2). This 
average energy is sensitive to the slope of the recoil spectrum, which is the quantity one 
needs to reconstruct fi{v). The statistical error of the reconstruction of fi{v) has been 
analyzed. Unfortunately it has been found that several hundred events will be needed 
for the method to be able to extract meaningful information on fi{v). This is partly 
due to the fact that fi{v) is normalized, i.e., only the shape of this distribution contains 
meaningful information, and partly because this shape depends on the slope of the recoil 
spectrum, which is intrinsically difficult to determine. 

Meanwhile, a method for determining the moments of fi{v) has also been presented in 
Subsec. 4.2.4. Numerical simulation shows that very rare events with large recoil energies 
contribute significantly more to the higher moments. Nevertheless, because this method 
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uses the whole experimental data together to determine the moments of /i (v) , it has been 
found that, based only on the first two or three moments, some non-trivial information 
can already be extracted from 0(20) events. 

As noted earlier, one needs to know the WIMP mass for the reconstruction of 
(the moments of) the velocity distribution. Moreover, although in well-motivated WIMP 
models can be determined with high accuracy from future collider data, we will have 
to check experimentally that the particles produced at colliders are in fact the same ones 
seen in direct Dark Matter detection experiments which form the Galactic halo. 

In Sec. 4.3 I have developed a method for (self-) determining based on the deter- 
mination of the moments of fi{v) by combining two (or more) experiments with different 
detector materials. The numerical analysis shows that, the larger the mass difference 
between two target nuclei, the smaller the statistical error will be. Hence, the combina- 
tions of two semiconductor detectors: Si and Ge and of two liquid noble gas detectors: 
Ar and Xe should be good choices. Meanwhile, due to the maximal measuring energy 
of the detector, there will be a deviation of the reproduced WIMP mass from the true 
one as WIMP masses > 60 GeV/c^. However, the numerical analysis shows also that, 
for WIMP masses < 100 GeV/c^ some meaningful information on the WIMP mass can 
already be extracted from C(50) total (each experiment C(25)) events. 

At the first step I have ignored the annual modulation of the WIMP flux. Given 
the large statistical errors expected in the foreseeable future, this is a reasonable first 
approximation. However, for the future detectors with strongly improved sensitivity and 
(very) large target mass (large exposure) , the formulae and methods have to be extended 
to allow for an annual modulation of the event rate. Hence, in Sec. 5.1 I have discussed 
the extension of the reconstruction of the velocity distribution by taking into account 
the time dependence of the recoil spectrum. The analysis shows that the two important 
observables for reconstructing the velocity distribution function: the event number and 
the average recoil energy measured in a given bin (or window), can be obtained as the 
annual average of the total event number (per day) and the average value of the recoil 
energies measured in experiments operated over some whole years. 

Moreover, in Subscc. 5.2.2 I have presented a method for reconstructing the ratio of 
the modulated amplitude of the velocity distribution to the time-averaged one. The only 
information which one needs is the measured recoil energies and their measuring times. 
This reconstruction is independent of the WIMP mass and can be done even if we can not 
obtain (enough) events in the high energy range. Hence, before the sensitivity of detectors 
can be improved to offer enough data until high energy range, reconstructing this ratio 
directly from experimental data and comparing it with the theoretical predictions might 
be the best possibility to test the different models of the halo Dark Matter. 

Furthermore, in Subsec. 5.2.1, I have given an alternative, and also better way to 
check whether the annual modulation of the event rate exists and thereby test models of 
the Dark Matter halo. The main advantage of this test is that, instead of (traditionally) 
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comparing the numbers of collected signal events in different, short time intervals in one 
year, one can now use information, i.e., the measured recoil energies and their measuring 
times, from all signal events collected in one or even several years together. For the case 
that the background events dominate the whole data set, this test might be still useful, 
if one expects that the background is (almost) time independent. 

The analyses of this work arc based on several simplifying assumptions. First, all 
experimental systematic uncertainties, as well as the uncertainty on the measurement 
of the recoil energy Q have been ignored. This should be a quite good approximation, 
given that we will have to live with quite large statistical uncertainties in the foreseeable 
future. Recall that, as shown in Sees. 3.7 to 3.9, not a single WIMP event has as yet 
been unambiguously recorded. 

I have also assumed that the detector consists of a single isotope. This is quite realistic 
for the current semiconductor (Si or Ge) detectors. On the other hand, for detectors 
containing more than one nucleus, by simultaneously measuring two signals, one might 
be able to tell on an event-by-event basis which kind of nucleus has been struck (see 
Subsec. 3.7.2). In this case, the methods can be applied straightforwardly to the separate 
sub-spectra. 

The analyses treat each recorded event as signal, i.e., background has been ignored 
altogether. At least after introducing a lower cut Qthre on the recoil energy, this may in 
fact not be unrealistic for modern detectors, which contain cosmic ray veto and neutron 
shielding systems (described in Subsec. 3.6). Background subtraction should be relatively 
straightforward when fitting some function to the data, which would allow to use the 
expressions given in Sec. 4.1. It should also be feasible in the method described in 
Sec. 4.2, if its effect on the average Q- values in the bins can be determined; in particular, 
an approximately flat (Q-independent) background would not change the slope of the 
recoil spectrum. 

In summary, a theoretical exploration of studying what direct Dark Matter detection 
experiments can teach us about the properties of Dark Matter particles in our Galac- 
tic neighborhood, e.g., their velocity distribution and their mass, the so-called "WIMP 
astronomy", has been started. However, the analyses show that this will require sub- 
stantial data samples. Hopefully this work will encourage our experimental colleagues to 
plan future experiments well beyond the stage of "merely" detecting Dark Matter. On 
the other hand, due to the significantly reduced condition (less than 100 events) for ex- 
tracting meaningful information on the WIMP mass by means of data from direct Dark 
Matter detection experiments, a championship for finding new particle(s) between the 
collaborations of direct Dark Matter detection and that of coUider experiments has also 
been started. 
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Appendix A 



Expression of the Velocity 
Distribution of WIMPs 



In this chapter I discuss at first some properties of the auxihary function Fi{v) defined 
in Eq.(4.1). Then I show two different approaches to find out the expression of fi{v) in 
Eq.(4.8), and derive the normahzation constant in Eq.(4.9) as well as the expression 
of moments of fi{v), (v^), in Eq.(4.10). 

A.l Properties of Fi{v) defined in Eq.(4.1) 

First, according to the definition of Fi{v) in Eq.(4.1) and noting that the velocity 
distribution function fi{v) can not be negative: 

/i(^)>0, (A.l) 



I have 

dF.iv) h{v) 



>0. (A.2) 



dv V 

This means that Fi{v) increases monotonically with v. Second, fi{v) must vanish as v 
approaches infinity: 

00)^0, (4.3) 

since WIMPs (as candidate for CDM) in today's Universe move quite slowly, then I have 



dFiiv) 



dv 



0. (4.4) 



v—^oo 



This means that Fi{v) must approach a finite constant Fi^oo as v — > 00. On the other 
hand, the three-dimensional velocity distribution function of WIMPs, f{v), must be 
bounded: 

f{v) ^IM^oo, (A.3) 

11^ 
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Figure A.l: Sketch of the auxiUary function Fi{v) defined in Eq.(4.1). 



where I have used 



J f{v)d\^ j f{v)v^dvdn^ J h{v)dv. 



Hence, 



dFAv) 



dv 



fiiv) 



v=0 



v=0 



0. 



(A.4) 



v=0 



This means that Fi{v) also approaches a constant Fi^o as w — > 0. Actually, Fi^o can be 
set to without loss of generality, since Eq.(4.1) defines Fi{v) only up to an additional 
constant. A sketch of the auxiliary function Fi{v) is given in Fig. A.l. 



A. 2 Derivations for fi{v) in Eq.(4.8) 

According to Eq.(3.10), I have 



dVr, 



a V-n 



dQ 2Q ' 



namely, 



dQ 2Q 



dVr, 



Differentiating both sides of Eq.(4.2) and using Eq.(4.1), one can obtain that 



/l(^^min) dFi{v^ 



d 



dvr, 



A (if n 



'dR' 
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^ min' 



(A.5) 



(A.6) 



~A\dQ 



dQ 



dv. 



1 



mm ^ 

dR\ 



^ min' 



J ) Q=t,2 . /q,2 
^ min' 



(4.5) 



namely, 



1 f 2Q 



1 



Q=i>2. loP- 

^ mm/ 

I (i / diC 



F{Q) \dQ J \dQ J dQ \dQ 



(A.7) 



min' 



On the other hand, according to Leibnitz's Rule for Differentiation of Integrals: 



d_ 
~dt 



bit) 

a(t) 



F{x, t) dx 



m 

a{t) 



dF{x,t) 
di 



dx + 



(A.8) 



one can also differentiate both sides in Eq.(3.12) with respect to Q directly and obtain 



d fdR 
dQ [dQ 

d 



—\af^{q) r 

dQ I -'Vmin 



h{v) 



'dF\Q) 


poo 


\fiiv) 


dQ 


•^^min 


V 



'dT 
dQ. 



2F{Q) 
2 ( dF\ fdR 



dv 



dv + AF^Q) 



' dVr, 



'dR" 



F\Q) \dQ^ 



-A 



F\Q) 



F{Q) \dQ) \dQ 
Therefore, it can be also found that 
1 [ 2Q 



2Q 



-AF\Q) 

fl ^^min ~ 



V 

'fi{n 



dQ , 

mmy / ^min 



2Q 



aJQ] . 



/l(^^min) 



_ 2 fdF\ fdR\ d (dR\ 

A \ FHQ) [ \dQ)\dQ)~'dQ\dQ)\\ . , 

^ min' 



(A.7) 



A. 3 Normalization constant and moments of fi{v) 

Since 

v = a^, (A.9) 
I have 



dv 



a 



dQ. 



(A.IO) 
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Prom Eq.(4.8) and according to the normalization condition in Eq.(4.7), I can get 

/ fi{v)dv 
Jo 

dQ 



1 dR' 
FHQ)\dQ^ 



a] 







= 1, 



Q 

oo 



^' dQ 

1 



'dR' 



F\Q) \dQ^ 



dR 



J 



^/Q 



F^Q) \dQ 
1 fdR\ 



dQ 

1 |.oo 1 

2J0 



dR' 



F\Q) \dQ^ 



dQ 



dQ 



where I have used the conditions: 

dR 

Q— »oo 



dQ 



0, 



and 



dR 
'dQ 



7^ 00 . 



(A.ll) 



(A.12) 



(A.13) 



Eq.(4.9) follows immediately from Eq.(A.ll). 

Using Eqs.(A.9), (A. 10) and integration by parts, I can also find the moments of fi{v), 
defined in Eq.(4.10) with a lower cut-off Qthre on the energy transfer, as follows: 



VminiQthie) 



v'^fi{v) dv 



/■OO 
JQth. 



aJQ] 1-2Q 



dQ 



'dR" 



a 



+1 

Qthre 

r n("+i)/2 

n+l ) V 



g(n+l)/2 . _^ 

dQ 



F\Q) \dQ^ 
d \ I ( dR 



F\Q) \dQ^ 



2VQ, 
dQ 



dQ 



ahre 



dR 

^^(Qthre) KdQ 



Q—Qthii 

n + l 



^ •'Qthre 

This reproduces Eqs.(4.11) and (4.12) in Sec. 4.1. 



'dR' 



[F^iQ) \dQj\ 



dQ 



(A.14) 
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Appendix B 



Moments of the Velocity 



Distribution of WIMPs 



In this chapter I derive at first the first and second moments of fi{v), i.e., the mean 
velocity and the velocity dispersion of WIMPs, for both of the two simplest semi-reahstic 
halo models discussed in Subsecs. 3.1.3 and 3.2.1. Then, as tests of the formulae for 
reconstructing fi{v) and determining I use the reduced spectra given in Eqs.(3.21') 
and (3.31') in Sec. 4.1 and the expressions for fi{v) in Eqs.(4.8) and (4.9) as well as for 
the moments of fi{v) in Eqs.(4.11) and (4.12) to obtain the same results. 



For a simple isothermal Maxwellian halo, the normalized one-dimensional velocity 
distribution function has been given as 



B.l Calculating (v) and {v'^) from fi{v) 



B.l.l From /i,Gau(^) given in Eq.(3.20) 




(3.20) 



One can find directly that 




(B.l) 
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where I have used 



fOO 

^0 



e-"^ dx 



Un+1) 



2a("+i)/2 ' 



and 



r(m + 1) = mr(m) , 



for m = 0, |, 1, |, 2, • • •. Hence, using 
r(n + l) 



and 



r(n + i) 



_ l-3---(2n-l)V7r 



2" 



for n = 0, 1, 2, 3, • • •, Eqs.(3.22) and (3.23) can be obtained directly. 

B.1.2 From /i,sh(^^) given in Eq.(3.29) 

When we take into account the orbital motion of the Solar system around the Galaxy, 
the velocity distribution function should be modified to 

1 / V 



First, I have 

roo 

{v)sh = / vfi,sh{v,Ve)dv 
J 



(3.29) 



1 



V ■ V 



1 



dv 



vq Jo vq Jo 



Define 



u± = 



^(Vi,_-Vi,+). 
^/nve^ ' 

V ± fp. 



^^0 



(B.2) 



(B.3a) 



I.e., 



V = VoU± ^ Ve , 



(B.3b) 
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it can be found that 
Vi,- 



and 



1 /■°° 
— / v'e 



vo Jo 
1 



— / {voU- + Ve) e (vodu-) 

Vq J -lie. 



u_e + 2veVo / M-e (iw 



Jve 



Vi,+ = 



1 r'^ 
Vq Jo 



Jo 

1 f'-'° 2 



u\e "+ du+ - 2veVo / w+e "+ + / e "+ du+ , 



Ve 



where I have defined 



^^0 



Combining Eqs.(B.4a) and (B.4b), one can get 
Vi,- - Vi,+ 



= 2^0^ (i;ee-^e) + erf (i;,) + 2^, 

= v^v,e-^'^/^l + ^ (I + ^e) (^) ' 
where I have used the definition of error function 

2 rx 



v^e-^- + 2v, 



TT 



^ I erf(^e) 



erf(x) = —= f e ^^dt, 

VTT Jo 



and 



Hence, the mean velocity of WIMPs for a shifted MaxweUian halo discussed 
sec. 3.2.1 can be found as 



(v) 



sh 



+ 



2v, 



+ Ve erf 
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Meanwhile, 



(^^^)sh = / v'^fi^si,{v,Ve)dv 
Jo 



1 



1 



v"^ ■ V 







dv 



Jo Vn Jo 



vo 



= ^(V2,--V2,+). 

Using Eqs.(B.3a), (B.3b), and (B.5), it can be found that 

1 f'-"-' 



(B.7) 



Vo_ = 



r v^e-^^'-^'^^'/^'o dv 
Jo 

/OO 2 
(voM- + Ve)^e""- {vodu-) 
-Ve 

/OO /" OO 

-■Ue J — Ve 

-Ve J —Ve 



Vq Jo 
1 /■'^ 

f 

(•OO 



/■OC 



u^e~^- du- + Sue^o 







+ 3VgVo / w_e du- + V, 

J tie 



(^''^ e""- + e""- du-^ ,(B.8a) 



and 



= 1 rv3g-(.+.e)V-g(^^ 

1 /""^ Q 2 

f -^J^e 

Jve J V, 



u±e "+ du^ 



(B.8b) 



Combining them, one can get 



Vo._ - V 



2,+ 



6VeVQ 



M^e-"' du + 2vl I e""^ du 



I 

Jo 



2v: 



^^fT^vl. 



(B.9) 



Therefore, the mean velocity of WIMPs for a shifted Maxwellian halo can be found as 

'3 



2 , „,2 
e 



Wo +t;, 



(3.33) 
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Finally, according to Leibnitz's Rule for Differentiation of Integrals given in Eq.(A.8), 
one has 



dx 



erf (x) 



2 
7^ 



d r 
dx 



^0 



dt 



(B.IO) 



Then it is easily to prove that, for Ve <^ Vq, i-e., Vf. <^ 1, Eqs.(3.32) and (3.33) will reduce 
to Eqs.(3.22) to (3.23). 

B.2 Calculating fi(v)^ (v), and {v'^) from dR/dQ 

B.2.1 From {dR/dQ)Ga^ in Eq.(3.21) 

Substituting Eq.(3.2r) in Sec. 4.1 into Eq.(4.8), I have 

'dR\ 11 

Gau- 



/l,Gau(^^) = A/bau { -2Q 



d 



Gau 



dQ 
d 



Gau 



(B.ll) 



Meanwhile, according to Eq.(4.9), the normalization constant Aoau can be found as 



Gau 



'I 



oo ][ 



'dR\ 



\F\Q)\dQ)^^^ 



a Jo 



dQ 



dq 



a \ a 2 ^ 
2 



Here, for simplicity, I have defined ^: 

Q = q\ 
and then 

dQ — 2qdq. 



(B.12) 



(B.13a) 



(B.13b) 



Substituting Eq.(B.12) into Eq.(B.ll), the normalized one-dimensional velocity distri- 
bution function /i,Gau('i') in Eq.(3.20) can be obtained directly. Moreover, substituting 



will use this definition in this section and the next chapter. Please do not confuse with the 
transferred 3- momentum in Eq.(3.7). 
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Eq.(3.21') into Eqs.(4.12) and (4.11) and using the normalization constant in Eq.(B.12), 
one can get 



a 



n+l ' 



{n+1 

,,n+l 



Jo 



1 fdR\ 



dQ 



a' 



M3au(n + 1) ( ^1 J^^ q--^e-^'iy-o {2qdq) 



'n + l' 



2\a 



|(n + l) 



|(n+l) 



TT 



Note that I have set Qthre — here. 



(B.l) 



B.2.2 From {dR/dQ)sh in Eq.(3.31) 

According to Eq.(B.lO), one can obtain that 



_d_ 
dQ 



erf( °^^^^- ) 



1 /a 



Then, substituting Eq.(3.31') in Sec. 4.1 into Eq.(4.8), I have 

/l,sh(^.^e) 



1 fdR" 
F^Q)[dQ^ 



shJ J Q=v'^/a'^ 



^ '^"^ ^ ^ -[{aV^+Ve)/vof _ -[{aVQ-Ve)/vo] 



=(-) 



sh 



I'D/ 



(B.14) 



(B.15) 



Meanwhile, as done for Acau in Eq.(B.12), one can use Eqs.(B.13a) and (B.13b) and find 
that 

-1 



A4 



1 /""^ 



sh 



a JO 



1 ( 



1 



dq 



a yjQ 
1 



eri — eri 



1^0 



^^0 



dq, 



a 



Vo(oo)-Vo(0) 



(B.16) 
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where I have defined 

Vo(g) ^ / 

Define 

aq ± V, 



aq + Ve\ , (aq - Ve 
eri — eri 



^^0 



^^0 



dq = Vo,+ (g) - 



s± = 



I.e., 



VoS± =F Ve 

a 



it can be found that 



erf(s-|- 

— — erf ds 
a J 



ds- 



Vo 

a 



s±erf(s±) + e '± 



H,±|)erf(.J.-L(;^,e- 



s2 



where I have used 



f 1 _ 2 

/ erf(a;) dx = a;erf(a;) H — ^ e ^ . 

J y/TT 

Substituting Vo,±(g) in Eq.(B.19) into Eq.(B.17), I can get 



Vo(?) = q erf(s+) — erf(s_) + — erf(s+) + erf (s_ 



+ 



Now note that, as g — * oo. 



Vo(g ^ oo) = ( - ) ' 



since 



and 



s±{q — > oo) — > oo , 



2 f°° 2 

erf(oo) = —;= / dt^l. 
■\/tt Jo 



While, since as = 0, 

s±(0) = ±^ = ±{;e, 

^"0 
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where I have used the definition in Eq.(B.18a), and 

erf(— x) = — erf(a;) , 
it can be found that 

Vo(0) = 0. 



(B.24) 



Substituting Eqs.(B.22) and (B.24) into Eq.(B.16), the normahzation constant AQh can 
be found as 

-I -1 



sh 



a 



1 

2vl 



(B.25) 



Then I can obtain the normahzed velocity distribution function in Eq.(3.29) directly. 
Meanwhile, substituting Eq.(3.31') into Eq.(4.11) (Qthre — 0) and using Eqs.(B.13a) and 
(B.13b), I have 



{v)sh = Kh ■ O? / 

1 ^ ("^ 
^ — .2o?\ q 
ZVf. Jo 



1 fdR\ 

aq + Ve 



erf 



^^0 



— erf 



{2qdq) 
aq — Ve 



Vq 



dq 



a 



Vi(oo)-Vi(0) 



where I have defined 



Vi 



erf 



aq + Ve 



erf 



aq — Ve 



dq = V^,+ {q)-V,,.{q). 



(B.26) 



(B.27) 



Using Eqs.(B.18a) and (B.18b), it can be found that 

Vi,±(g) 

'aq ± Ve 



^/gerf( 



a 



Vo 

VoS± T Ve 

a 



dq 

erf (s±) ds± 



= f — ) / s ±eif {s±) ds± T / erf (s±) (is 



VeVo 



a^ 



s±erf (s±) + — e *± 



1 /^^^2 

2\a 
1 fvo^^ 



2Ve 



s±erf H — — e 



TT 



-erf(s±) 



(B.28) 



where I have used 



J xeri{x) dx = - 



X 



^ ) erf(x) + xe 



TT 
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and 



VoS± ^ 2ve = {aq ± Vg) ^ 2ve = aq^Ve = VqS^ . 
Hence, I can get 



erf(s+) — erf(s_) 



From Eq.(B.22), it can be found easily that 

Vi(g^oo) = 0, 
and, from Eq.(B.23), 



Vi(0) 



a J 



(B.29) 



(B.30) 



(B.31) 



(B.32) 



Therefore, substituting these results into Eq.(B.26), I can obtain that 



a 



a J 



Similarly, 



sh 



1 

2v, 



2v, 



3 
2 

■3a^ 



2 /„2 



TT 



(3.32) 



a 



1 

erf 



1 (dR\ 



f aq + 



erf 



(2gdg) 
f aq — Ve 



dq 



- 3 K 
where I have defined 



V2(oo)-V2(0) 



V2(g) = / 



erf 



-"0 



— erf 



otq - Ve 



dg^V2,+ (g)-V2,-(9) 



(B.33) 



(B.34) 



Using Eqs.(B.18a) and (B.18b), it can be found that 

'Q:q ± Ve 



V2,±(g) = J q'eri 



dq 



= (5) / (""o^i ^ ^VeVoSi + erf (s±) 

y s^erf (s±) ds± =F y s±erf (s±) ds± + v^ j erf (s±) ds± 
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4erf(s±) + ^(4 + l)e-^± 



\ ) erf(s±) + ^s±e-"± 



I — 

1 _ 2 ' 
s±erf (s±) + — e *± 

TT 



1 / a \ 



erf(s±) 



where I have used 



x^erf(a:) dx = - x^Qxi{x) + ^(x^ + l) 
3 v"^ 



and 



'± ~ 2 ~2 



1 / T;og± ^ 

3 V Wo 
1 

3 Vvo 



3 



Hence, I can get 



V2(g) 



+ 



erf — erf 



+ — + — 

I 3 2 



erf (s+) 



+ 



Ve s+e 



+ Ve S-e~ 



e 4 - e-- 



From Eq.(B.22), it can be found easily that 



V2(g^00)=(^ 

a 



Ve_ 
3 



a 



and, from Eq.(B.22), 
V2(0) = 0. 

Therefore, substituting these results into Eq.(B.33), I can obtain that 



3 fa^' 

2 I Ve 



a"" 
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Appendix C 

Differential and Total Event Rates 



In this chapter I derive the differential and total event rates for the simple and shifted 
isothermal Maxwellian halo models from their velocity distribution functions given in 
Eqs.(3.20) and (3.29). The case for F'^{Q) fti 1 and the case with the exponential form 
factor F^^{Q) given in Eq.(3.15) will be considered. 

C.l Setting F\Q) ^ 1 

C.1.1 Starting with /i,Gau(^) given in Eq.(3.20) 

For a simple isothermal Maxwellian halo, the normalized one-dimensional velocity 
distribution function has been given as 



/l,Gau(^') 

I can get directly that 

/ 

Jv, 



(3.20) 



/l,Gau(t') 


dv — 




V 







dv 



2 



3 / 



ve 



min 

OO 



Using this result and Eqs.(3.10) and (3.12), one can obtain [dRj dQ;)Q^^ in Eq.(3.21) and 
then i?Gau(Qthre) in Eq.(3.24) easily when F'^{Q) has been neglected. 

C.1.2 Starting with /i,sh(^) given in Eq.(3.29) 

When we take into account the orbital motion of the Solar system around the Galaxy, 
the velocity distribution function should be modified to 



/l,sh(^^,^^e) 



1 



(3.29) 
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First, using Eqs.(B.3a) and (B.3b), it can be found that 

1 



poo 


/l,sh(^^) 


' ^min 


V 



dv 



1 



•^^"1X1111 



1 



P oo /* oo 



1 

2^ 



'^^+,min 

erfc(ti_,inin) - erfc(M+ 



erf (w+^min) - erf (m_ 



(C.2) 



where I have used the definition 

2 ,2 

erfc(a;) = —j= I e dt = 1 — erf (x) , 

y/TT Jx 

and 



U±,r> 



^min i 
^^0 



(C.3) 



Combining Eqs.(C.2) and (C.3) with Eqs.(3.10) and (3.12), one can obtain {dR/dQ)sh in 
Eq.(3.31). Moreover, by using Eqs.(B.13a) and (B.13b), one can find that 



Qthre 
OO 



\ vo J V ^'o 



dQ 



roc 

"'9th: 



aq + Ve\ faq-Ve 
erf — erf 



= 2 



Vo y \ Vo 

Vi(g ^ oo) - Vi(gfthre) 
2 r /I 



{2qdq) 



3 115-^^^ 



erf (5+) - erf (5_ 



Here I have defined 



(C.4) 



(C.5) 



and used Eqs.(B.27), (B.31), (B.30), (B.18a), and (3.35). Hence, for F^{Q) ^ 1, one can 
get -Rsh(Qthre) in Eq.(3.34) directly. 

C.2 Using Fl{Q) given in Eq.(3.15) 

C.2.1 Starting with {dR/dQ)Gan given in Eq.(3.21) 

Substituting the exponential form factor F^^{Q) given in Eq.(3.15) into Eq.(3.21), one 
can get 



'dR\ 
dQj 



Gau.ex 



yv^voj w^voj 
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where (5 has been defined in Eq.(3.28). Then it is easy to find that 

it:Gau,ex(gthre) = A H e'^'^/'^o^^ dQ 



(3.26) 



C.2.2 Starting with {dR/dQ)sh given in Eq.(3.31) 

Substituting the exponential form factor F^^{Q) into Eq.(3.31), one can get 



'dR\ 



-Q/Qo 



erf( "^+^- ) - erf( "^-^- ) 



(C.7) 



Consider 



J e-^/^° erf(^^^^) dQ 
^ -Qoe-'^/^^erff^^Si^l + f-) / ^ e-«/«°-[("^±-«)/''°]'dg , 

where I have used integration by parts and Eqs.(B.14). Using Eqs.(B.13a) and 
the integral of the second term on the right-hand side above can be found as 



(C.8) 

(B.13b), 



dq 



= V^(^)e-M^Kerf(^±/3i;.), (C.9) 



where 1 have used 



j ^-(„.^+6x+c) dx^\^ e^"'!'^-^) erf (^^x + ^) , 
and the definition in Eq.(B.5). Combining Eqs.(C.7) to (C.9), one can get 



-Rsh,cx (Qthrc) 



POO-Q / % 



m-^rrif^ \2veJ \1 — (3"^ ) 
X /e-(i-/3')"'Qth.eK'/3' 



erf (5+) - erf (5_) 



erf(T+) -erf(T_) 



where I have used the definitions in Eqs.(3.35) and (3.38) as well as 



<5o = 



vl( 

q;2 I 1-/32 



(3.37) 



(3.28') 
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Appendix D 

Some Old Attempts 



In this chapter I present some old attempts for reconstructing the velocity distribution 
function, eventually also for determining its moments. I describe also their disadvantages 
and problems. However, these unsuccessful attempts could perhaps inspire some new 
ideas. 



D.l Binning the data set 

The usually used choice for binning a data set is that every bin has that same width: 

Quisx Qmin ^j-^ 



bn = b^ 
and thus 

Qn — Q 



B 



n ]h . 

2/ 



(D.2) 



However, as discussed in Subsec. 4.2.2, using bins with linearly increasing widths can 
make the errors roughly equal: 



here the increment 5 satisfies 

5 — _ (<5max — Qmin — . 

Hence, for the n-th Q-bin, one has 

^ / .XL r(n-l)(n-2) 

and 

Qn,Taax. — Qmin ~l~ Tlh\ -\- 

This means that 



n(n — 1) 



8. 



Qn = Qmin + ( " 2 ) ^1 + 



fn-ll 



8. 



(4.33) 



(4.35) 



(D.3a) 



(D.3b) 



(4.34) 
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Moreover, one other choice for binning the data set is 

bn = hS^-' . (D.4) 
It is more comfortable if we choose 5 as the input parameter and then determine 61 as 

(D.5) 



(5 — 1 \ 
^ j (^Qmax Qmin 

Hence, for the n-th Q-bin, one has 



Qn.min Qimn 



6-1 



and 



" - 1 

Qn.max — Qmin ~l~ | ~J ~ ) ^1 



6-1 



This means that 



Qn — Qmin ~l~ 



5" + 



2(5-1) 



bi. 



(D.6a) 



(D.6b) 



(D.7) 



D.2 Reconstructing fi{v) without derivatives 

According to the expression of the differential event rate in Eq.(3.12), I have 



'dR\ 



AF' 



'iQn) / 



dv , 



where, from Eq.(3.10), 

Vfi Qn • 

Then it can be found that 

fVn+l 



(D.8) 



(D.9) 



j-Vn 
Jvn 



\fi{v)] 


dv = — 


V 


A 



1 fdR\ 



1 fdR\ 



F\Qn) \dQJo=o P\Qn+i) \dQ) 



Q=Qn+l, 



A„.. 



The mean value theorem of calculus implies 



fi{v) 



[D.IO) 



[D.ll) 



where Vn < Vn < Vn+i- Hence, I can let 

Vn = 0£nVn+l + (1 " Oin)Vn = t'n + 0£n{Vn+l " Vn) , 

and rewrite Eq.(D.ll) to 



Vn+1 - Vn^ 



1 



{Vn+l/Vn) - 1 
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(0 < < 1), (D.12) 



A„, (0<q;„<1). (D.13) 



Therefore, the error of fi{vn) can be given as 

1 



<7(A„) . 



_{Vn+l/Vn) - 1 

Usually, one sets an — \ and then it can be reduced to 

(T{h{Vn,an = 1/2)) = Ct(/i(v„+i/2)) = ^ [ !^!1±1^^ 

Here, from Eq.(D.lO), 



dQ)Q=Q„ 



1 



1 [ 1 /iv, 



+ 



(Nn+l 



a 



(IR 



1/2 



n+1 



1/2 



(D.14) 



(D.15) 



(D.16) 



where I have used the standard estimator for for dR/dQ at the point Q = (5„ in Eq.(4.15) 
and then its statistical error in Eq.(4.16). 

This method is straightforward. However, neither fi{vn) nor its statistical error is 
independent of the unknown constant A. Moreover, this method has an anti- correlation 
problem: An upward fluctuation of the counting rate in the n-th Q-bin will lead to too 
small /i in the n — 1-st v-hin, but tends to give too large /i in the nth v-hin. 



D.3 Average logarithmic slope 

As shown in Fig. 4.1, the theoretically predicted recoil spectrum is approximately 
exponential. And, as discussed in Subsec. 4.2.1, an exponential approximation can ap- 
proximate the recoil spectrum for a wider bin. Hence, 1 have considered the exponential 
ansatz in Eq.(4.18), but at beginning only naively combined with the standard estimator 
for dR/dQ at the point Q = Qn a.s 

= r^e^"^^-'^") , (D.17) 

where r„ = Nn/hn is the standard estimator given in Eq.(4.15). Define the slope of the 
straight line with two endpoints {QnM{dRl dQ)Q=Qu) and {Qn+iM{dRl dQ)Q=Qn+^ ^ 

_lnr„+i-lnr„ m 1Q^ 

kn,n-\-i = —^ Tn — ' ^ = 1,2, ■■■,B-1. (D.18) 

Then 1 can define an average slope for the function ln{dR/dQ)Qc^Q^ at the point Q = Qn, 
n = 2, 3, • • • , 5 - 1, as (see Fig. D.l) 

^ ^ kn-i,n + kn,n+i ^ 1 / lu r„ - lur^-i _^ lnr„+i - lnr„ \ ^^^^ 

2 2 y Qn — Qn-l Qn+l ~ Qn J 
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Figure D.l: Sketch of the average slope for the function h\{dR/ dQ^Qn^q^ sX Q = Qn, 
^n,ave, defined in Eq.(D.19a). 



but, at the point Q = Qi, I have defined 
_ In r2 — In ri 

"'l.ave = '^1,2 



(D.19b) 



Q2-Q1 

The statistical errors on /Cn^ave can be obtained directly from Eqs.(D.19a) and (D.19b) as 



^Qn Qn—1 Qn+1 Qn J -^n 

\ 2 



+ 



for n = 2, 3, ■ ■ , B — 1, and 



1 



+ 



1 



Qn ~ Qn-1 J ^n-l \Qn+l ~ Qn J ^n+l 



where I have used Eqs.(4.15) and (4.16) to get 



„2 
n 



(D.20a) 



(D.20b) 



(D.21) 



Moreover, the kn,ave given in Eqs.(D.19a) and (D.19b) are correlated. Hence, for /cn.ave in 
Eq.(D.19a), I have 
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COV (Ajjjjave; ^n+l,ave) 



1 



Qn Qn—1 Qn+1 Qn J \Qn+l Qn / 

1 \ / 1 1 



+ 



.Qn+l ~ Qn J \Qn+l ~ Qn Qn+2 ~ Qn+1 J ^n+l 



and 



COV {kji^avej ^n+2,ave) — 



4 \Qn+l — Qn) \Qn+2 — Qn+l) ^n+l ' 

while, for /ci,ave in Eq.(D.19b), I have 



(D.22a) 



(D.22b) 



COV (/Ci^ave) ^2,ave) — 



1 



1 

-T7- + 



1 



and 



GOV (/Ci,ave, h,s 



1 



^Q2-Qi) Q2-Qi\Q2-Qi Q3-Q2) N2 

(D.22c) 



(D.22d) 



2 \Q2-Q1) \Q3-Q2) N2 

Now I can begin to reconstruct the recoil spectrum. The basic idea is that I approx- 
imate the function \ia{dR/dQ) in each bin by a straight line lnr„,ave(<5) which has the 
slope kn,ave and passes through the point ((5„,lnr„) (see Fig. D.2): 

lnr„,avc(Q) - lnr„ 



Q-Qn 

Hence, I have 

in the n-th Q-bin: 

Qn-l+Qn — <" n <" D — Qn+Qn+1 

2 — Vn— _ _ ^;n+ — 2 ' 

and 

/ ^ d\ 

ri,ave(Q)=rie'=i--«-«i), 



'dR\ 

^d'QjQ^Q, 



in the first Q-bin: 

Qtl.re = Ql-<Q<Ql+=^^, 



(D.23) 



n = 2, 3, (D.24a) 



(D.25a) 



(D.24b) 



(D.25b) 
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Figure D.2: Sketch of the reconstructed segment of the function ln{dR/dQ) between 
{Qn-i + Qn)/^ and {Qn + Qn+l)/'2, In 

' n,ave iQ)- 



where Qthre is the threshold energy and the /cn,ave, = I, 2, ■ ■ ■ , B — 1, are given in 
Eqs.(D.19a) and (D.19b), respectively. Then, similar to Eq.(4.46), the velocity distribu- 
tion function fi{v) given in Eq.(4.8) can be reconstructed as 



/l,ave(^'n) = Mi^ 



(D.26) 



for n = 1, 2, ■ ■ ■ , B — 1. Here Vn is given in Eq.(D.9). 

The first problem with the expression in Eq.(D.26) is estimating the normalization 
constant A/lve- One possibility is inserting r„^ave(Q) given in Eqs.(D.24a) and (D.24b) 
into Eq.(4.9) directly: 

^ -1 



Qi- \fQ 



dQ 



(D.27) 



where Qi± are given in Eqs.(D.25a) and (D.25b). Similarly, defined in Eq.(4.12) and 
{dR/ dQ)Q=Q^^^^ in Eq.(4.11) can also be estimated as 



rQi+ 



dQ, 



and 



'dR^ 
dQ, 



ri,ave(<9thre)=rie'=i'-«(«*''--«^). 



(D.28) 



(D.29) 



ave,Q=Qthre 
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Q 



1 



n 



Qn + Qn+± 

2 



Figure D.3: Sketch of the elevation from {dR/ dQ)^e^\ Q=Q^ to r„ and from r„ to r* due 
to the concavity of the recoil curve, {dR/dQ)^ea.i, and that of the reconstructed spectrum, 
rn,iwe{Q), between (Qn-i + Qn)/2 and {Qn + Qn+i)/2, respectively. 

However, it should be pretty complicated to estimate the statistical errors of /i,ave('i^n) 
given in Eq.(D.26) with A/lve estimated in Eq.(D.27). 

The other serious problem with the ansatz in Eqs.(D.24a) and (D.24b) is that one 
must also consider a systematic error caused by using the exponential ansatz with the 
standard estimator r„ = Nn/h^. Suppose that {dR/dQ)j-ea.i is the real recoil spectrum and 
passes through the point (Qn, {dR/dQ)rea.\,Q=Q„) (see Fig. D.3). During the experiment 
we measure deposited energies and count the event rate, which is proportional to the 
area under the real recoil spectrum (see Eq.(4.20)), in the n-th Q-hin, and then estimate 
r„ = Nn/bn- However, because the recoil spectrum is concave, the estimator Trj^ IS a 
little larger than the real value, {dR/dQ)^ea.i,Q=Q„ (see r„ given in Eq.(4.22)). Define this 
elevation from {dR/dQ)rea.i,Q=Q„ to r„ as 



On the other hand, it is plausible to suppose that the reconstructed recoil spectrum 
rn,scve{Q) in Eq.(D.20a) and (D.20b) are approximately parallel to the real one {dR/dQ)rea.h 
thus I can estimate the elevation by 




n ' n,: 



n= 1, 2, B -1. (D.30) 



A, 



■r,n 




— r. 



n = 



1, 2, ■■■,5-1. 



(D.31) 
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Here r* a,ve can be calculated from r„,ave(Qn) as 

2 /■{Qn+Qn+l)/2 



n,ave ^ _ ^ 

2 



- / rn,a.ve{Q)dQ 

-1 "'(Qn-l+Qn)/2 



n-1 J{Qn-l 



(Qn+Qn+l)/2 



+Qn)/2 



(Q-Qn) 



_^n,ave(gn+l Qn—l)_ 

for n = 2, 3, • • • , S — 1, and for n — 1, 

2 /■(Qi+Q2)/2 



, / Qn+l-Qn ^ , ( Qn-Qn-\ \ 



, (D.32a) 



l,ave 



Ql + Q2 — 2(5thre JQ thre 

2 /■(Qi+Q2)/2 
2 

ri 



ri e 



fcl,ave(Q-Ql) 



dQ 



{—^ — -) _ g-fcl,ave(Ql-Qthre) 



(D.32b) 



.^l,ave(gi + Q2 — 2(5thre). 

Combining Eqs.(D.30) to (D.32b), the real value of the recoil spectrum {dR/dQ)reai at 
the point Q — Qn, n — 2, 3, ■ ■ • , B — 1, can be obtained (approximately) as 



^ n.ave 



= 2r„ 1 



1 



^n,ave(Qn+l Qn—1, 

, f Q„+i-Qn \ 
'^n.ave 1 9 I 



— e 



(D.33a) 



and at the point Q = Qi, 
ri,real ~ 2ri - r* 
= 2rJl 



^l,ave(Ql + Q2 — 2(5thre). 



X 



,ave 

^Q2_Qi) _ g-fci ,ave(Ql— Qthre) 



(D.33b) 



Note that the correction of r„ here is essentially the same as the expression of r„ in 
Eq.(4.22). 

Now one can replace r„ in Eqs.(D.19a) and (D.19b) by r„ real estimated by Eqs.(D.33a) 
and (D.33b) to get A;„ avcj and then substitute r„ i.cai and kn,sve into Eqs.(D.26) to (D.29) to 
reconstruct /i,ave('i^n) and so on. However, due to the dependence of r„ real on r„ and kn,a.ve, 
it is very complicated to modify even the statistical errors in Eqs.(D.20a) to (D.22d)! 

Moreover, this "average-logarithmic-slope" method has also the same "anti-correlation" 
problem as the method described in the previous section. An upward fluctuation in the 
n-th Q-bin leads to a too small slope and a too large slope kn^n+i, even though 

the fluctuation of the "average" slope kn^avc = (^n-i,n + kn,n+i)/2 could be more or less 
decreased and the value of /cn,ave should be not very bad. 
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Furthermore, as shown in e.g., Eqs.(D.18), (D.24a), and (D.30), from B bins one 
can get S - 1 A;„,ave, '^n,ave(Q), and r„,reai- Then, after one replaces r„ in Eq.(D.18) by 
rn,Tea.i and runs the whole process from Eq.(D.18) to Eq.(D.29), one can have only B — 2 
fi,ave{vn) givcu by Eq.(D.26). However, as shown in Figs. 4.2, with 500 (or even 5000) 
events, one has only 4 (or 8) bins to use. ^ Hence, it can not be allowed to lose 2 bins 
(points) more! 



D.4 Linear approximations of {dQ/dlVjQ^Q^ 

As noted in the beginning of Sec. 4.2, according to the expression in Eq.(4.8), one 
needs not only an estimator for dR/dQ aX Q = Qn but also one for the slope of the recoil 
spectrum to reconstruct the velocity distribution. A rather crude estimator of this slope 
is 



/ fdR" 
dQ VdQ^ 



Q=Qn 



(D.34) 



where Nn^q^q^ and Nn,q<q„ are the numbers of events in bin n which have measured 

recoil energy Q larger and smaller than Q„, respectively. This estimator is rather crude, 
since it only uses the information in which half of its bin a given event falls. 

It is clear intuitively that an estimator that makes use of the exact Q-value of each 
event should be better. This can e.g., be obtained from the average Q- value in a given 
bin. Taylor-expanding dR/dQ around Q = Qn, keeping terms up to linear order, gives 



'dR' 
dQ, 



d /dR" 
dQ [dQ, 



(Q Qn)^n 



Q=Qn 



(D.35) 



(D.36) 



Using this linear approximation for the recoil spectrum, one can find 

rQn+b-a/2 |- , 
Nn^ Tn+iQ - Qn)Sn dQ = r„6n , 

(of course, this reproduces the standard estimator in Eq.(4.15)) and the average value of 
the recoil energies in the n-th Q-hin: 



(•Qn+bn/'i 



1 /■<9n+ . 

Q - Qn\n = TT / {Q - Qn) rn + {Q - Qn)Sn dQ 

J^n JQn-bn/2 L ^ 

Hence, an improved estimator of the slope of dR/dQ at Q — Qn is 



(D.37) 



S2,n = 



12r,Q - Qn 



(D.38) 



^Not 5 or 10 bins, because the last one or two bin are almost empty, see the discussion in Subsec. 4.2.3. 
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According to the definition of si,n in Eq.(D.34) it can be found that 



(7^(si,n) 



1 


2 


{hnl2)\ 





16 



(at, 



hi 



n,Q>Q„ + Nn,Q<Q, 



(D.39) 



where I have used Eqs.(4.17) and (4.15). On the other hand, according to the expression 
of S2,n in Eq.(D.38), and treating the number of events and the average Q- value in a given 
bin as two independent variables, one can obtain that 

2 /-, „ \ 2 



2/_ \ ( 12Q Qn|n 



O- (S2,n) = 



62 



Here I have used Eq.(4.16), the definition in Eq.(4.24) with the linear approximation in 
Eq.(D.35), 



{Q-Qn)\ = Tr 



1 rQr.+b-a/'i 



{Q - Qnf Tn+iQ- Qn)Sn dQ = ^ , 



12 



and 



C (^Q Qn\nj — {Q Qn)'^\n Q Qn\n 



(D.41) 



(D.42) 



This simple calculation shows that the estimator S2,n given in Eq.(D.38) indeed has a 
smaller statistical error than the crude estimator in Eq.(D.34) by a factor of ^J^J^. 



D.5 Using the exponential ansatz in Eq.(4.18) 

In App. D.3 I have used an exponential approximation with the standard estimator 
r„ to reconstruct the recoil spectrum. A correction due to the approximately exponential 
form of the recoil spectrum has also been discussed. The use of the average logarithmic 
slope kn,a-ve Combined with the correction is very complicated, especially for the error 
analysis. However, it is clear that an exponential approximation can approximate the 
recoil spectrum much better than a linear one. On the other hand, in the previous section 
I have introduced the use of the exact Q- value of each event. The analysis done with 
two linear approximations has shown that the statistical error can be strongly reduced. 
Hence, it is pretty straightforwardly to combine these two techniques and their advantages 
together. 

By using an exponential ansatz for the recoil spectrum in each Q-bin, combining with 
a prefactor which can be adjusted by the event number in this bin, and then estimating 
the logarithmic slope by the average value of the recoil energies measured in this bin, 

^Strictly speaking, the denominator should be Nn — 1 as I used in Eq.(4.32). 
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one can already obtain the expressions given in Subsec. 4.2.1. Substituting the first 
expression of the exponential ansatz in Eq.(4.18) into Eq.(4.8) with the logarithmic slope 
kn estimated by Eq.(4.23), the velocity distribution function can be reconstructed as 



fl,Qi'"n) = ^fQ 



kn 



Q=Qn 



;d.43) 



where f„ is given in Eq.(D.9). This expression is already almost the same as the expression 
given in Eq.(4.46), except the central point Qn in the n-th Q-bin has been used here 
instead of the shifted point Qs^^ in the /x-th Q-window, and thus 1 have used fVi instead 
of here. Moreover, I have to determine the normalization constant Mq here. It can be 
done by Eq.(D.27) with replacing rj_avc(Q) by f^e'^"^'^'"'^"): 



a 



where Qn± have been given as 

Qn ~2' ~ Qn— — Q — Qn+ — Qn ~l~ 



(D.44) 



(D.45) 



However, in order to estimate the statistical error of fiqivn) in Eq.(D.43) more easily, I 
have defined 



fl,nd'"^ = 

and 

/l,n,Q(^) = 2 

in the n-th t'-bin: 



F\Q) 



F\Q) 



dlnF^iQ) 
dQ 



kn 



d\nF\Q) 
dQ 



kn 



/l,n,Q(^) 



Q 



aVQn- = Vn-<V < Vn+ = O^^JQn^ , 



;d.46) 



(D.47) 



(D.48) 



for n = 1, 2, 3, • • • , B. The normalization condition in Eq.(4.7) can be rewritten as 



1=1 Vi = l ■'^i- 

Then the normalization constant can be obtained directly by 



(D.49) 



B 



i=l 



E / kid"") 



(D.50) 



and find^^ Eq.(D.46) can be rewritten as 



/l,n,Q(^) 



.1=1 •''"i- 



/l,n,Q(^) 



/l,n,Q(^) 



(D.51) 



131 



Here I have defined 

B 



i = l ■'''i- 

Moreover, it is reasonable to define the n-th moment of fingiv) as 



:D.52) 



(D.53) 



Furthermore, in order to estimate the statistical errors of fi^qi'^) ^^^1 {v"')q given in 
Eqs.(D.51) and (D.53), I have denoted first the independent variables of /i „ ^('y) defined 
in Eq.(D.47) as x^j, where the subscript u stands for different species of variable, and 
j = 1, 2, 3, ■ ■ ■ , B stands for the bins. Meanwhile, I have assumed that, in the j- 
th Q-bin (and then also in the j-th v-hva), the error of each of these variables x^^j is 
approximately equal. Hence, I can use its value at the point Q — Qj, defined as 5xi,,j, for 
the whole j-th Q- or v-hin. Prom the expression of /i„q(v) in Eq.(D.51), its statistical 
error can be found directly as 



1/ j=l 
1 



dx 



C4 /— 
^0 1/ 



a'ix^j) . (D.54) 



Here I have defined 



.f 



and 



'S'l/JlA = E 



i = l •' "i 



dv . 



;d.55) 



(D.56) 



Similarly, from the definition of {v"')-q in Eq.(D.53), it can be found that 



u j=l 



dx^.j 



a {x^,j) 



I B 2 

<^ E E ('^0 ^J^,3\n — Sn Su,j;o) {x^^j) 
^0 V j=x 



(D.57) 



According to the expression of jx^n^{v) in Eq.(D.47) and and A;„ given in Eqs.(4.22) 
and (4.23), the independent variables of jx^qiv) should be chosen as 



(D.58) 



132 



with a'^{Nn) and a'^{kn) given in Eqs.(4.17) and (4.29). Since finqi^) depends only on 
Nn and km Eqs.(D.54) and (D.57) can be reduced to 



1 r ~ ~ 



and 



a'{x,j) , (D.59) 



(D.60) 



Here I have defined 



(D.61) 



dlnF^(Q) 



H -1 



and, from Eq.(D.56), 



dv . 



(D.62) 



(D.63) 



D.6 Introducing the average value of Q^lF'^{Q) 

The method presented in the previous section has two disadvantages. First, the 
estimator of Nq in Eq.(D.50) is the sum of several integrals, this makes the estimation 
complicated. Second, and also the worse disadvantage, by using {v^)-q given in Eq.(D.53) 
with 5";^ defined in Eq.(D.52), one has to know /i „ defined in Eq.(D.47), i.e., /i „ ^('y) 
defined in Eq.(D.46). It is not only comphcated but also loses the advantage of the 
expressions in Eqs.(4.11) and (4.12), by which one can evaluate the moments of fi{v) 
without knowing the functional form of fi{v). This problem comes essentially from the 
estimator of N'q in Eq.(D.50) obtained from the normalization condition in Eq.(4.7). 
Hence, one needs a new estimator for the normalization constant. 

Similar to the use of the moments of the recoil spectrum in Eqs.(4.23) and (4.24), I 
have defined an average value of Q^/F'^{Q) for all events in the n-th Q-bin: 



-/ 

N„. Jc 



Qn+b„/2 Q> 



dR 



:^HQ = ^E 



Then, for all recorded events in the sample, I can use 



= 8. 



(D.64) 



(D.65) 
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Note that the recoil spectrum dR/dQ here is not specified to the exponential ansatz 
{dR/dQ)n in Eq.(4.18). Using the definition in Eq.(D.65), N in Eq.(4.9) and 4 in 
Eq.(4.12) can be estimated by 

-1 

(D.66) 



a \ 62, 



-i/2,tot/ a J 



and 



2,(n-l)/2,tot 



E ^nS2,(n-l)/2,n ■ 



;d.67) 



n=l 



The expressions in Eqs.(D.66) and (D.67) are already essentially the same as the expres- 
sions given in Eqs.(4.48) and (4.54), respectively. 

Now replacing jV'g in Eq.(D.43) by given in Eq.(D.66), the reconstructed velocity 
distribution function at point v — Vn, /i,n,Q('^n), in Eq.(D.43) can be expressed simply as 



Q=Qn 



f 



l,n 



(D.68) 



tot , 



Here, similar to /i„q('1') defined in Eq.(D.47), I have defined 



fl,n — 



_d_ 

dQ 



In F\Q) 



Q — Qn 



[Dm) 



Moreover, the n-th moment of the velocity distribution function, (v"), determined by 
Eq.(4.11) can now be expressed as (with Qthre = 0) 



(^;") = a"(n + 1) 



'5'2,(n-l)/2,tot 
, 'S'2-l/2,tot , 



(D.70) 



By means of this expression, one can finally estimate the moments of the velocity distri- 
bution function directly from the experimental data given in Eq.(4.14) without knowing 
the exact form of /i(fn). Actually, according to Eq.(D.67), the expression of {v"') given in 
Eq.(D.70) is exactly the same as that given in Eq.(4.66), or the general form in Eq.(4.52) 
with Qthrc = 0. 

On the other hand, before beginning to calculate the statistical errors of /i(fn) and 
(f") in Eqs.(D.68) and (D.70), one must pay some special attention with the variables 
involved in their expressions. According to the expression of in Eq.(D.69) and of f„ 
and kn in Eqs.(4.22) and (4.23), is a function of only two variables: Nn and Q — Qn\n'i 
while, according to the definition in Eq.(D.67), 5'2,A,tot is a function of 2B variables: Nn 
and S2,\,n for all n = 1, 2, ■ ■ ■ ,B. Hence, /i(fn) and (w") depend on 2B + 2 and 3B 
variables, respectively. 

Then, from Eq.(D.68), it can be found that 



_£fiin 
dQ 



=^ = -fl{v„) (Q -Qn\n + Kn 



(D.71a) 
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where Kn has been defined in Eq.(5.31) and I have used dQ — Qn\n/dkn in Eq.(4.30); for 
m = 1, 2, • ■ ■ , B, one has 



-1/2, m 



= fl{Vn) 



Nm S2, -l/2,tot. 



and 



= -/iK) 



dSo^ 



2,-l/2,m 



'2,-l/2,tot, 



(D.71b) 



(D.71c) 



Then the statistical error of fi{vn) estimated in Eq.(D.68) can be expressed as 



Q Qn\n + -^n) 



Cx'(g-Qnln) + 



Nr,. 



+ 2 



5'(«n) 

'S'4,-l,tot 25'2,-l/2,n 



'S'2,l/2,(i ^ Q/i'-92,-l/2,n 



5. 



2,-l/2,tot 



'S'|-l/2,tot 'S'2-l/2,tot 



(D.72) 



Here I have used Eq.(4.17) for cr^(A^„), and, for simphcity, set ^ 1 for all m in order 
to use 



and 



COV (Q — Qn\n, S2,X,n) — ('S'2,A+l,n " QnS2,X,n) , 
COV (S'2,A,n, S2,p,n) = ('S'4,A+p,n " 'S'2,A,n S2,p,n) , 



with the definition 



'4,A,n 



1 -Afn /OA 



(D.73) 



(D.74) 



(D.75) 



and then 



'S'4,A,tot = ^nS4:,X,n i 



(D.76) 



n=l 



see Eqs.(D.64) and (D.65). Meanwhile, from the expression of (i'") given in Eq.(D.70), 
it can be found that 



d{v") 



dS2 ,(n— l)/2,m 



dS. 



2,-l/2,m 



,'S'2,(n-l)/2,tot, 

^ m 1 



2,-l/2,tot, 



(D.77a) 
(D.77b) 



135 



and 



S2,{n-l)/2,m 
,5'2,(n-l)/2,tot 



52,-l/2,m 
'2,-l/2,tot, 



fD.77c) 



Hence, the statistical error of (w") can be expressed as 



S4,n-l,tot _|_ 'S'4^_i^tot 



2S. 



4,(n-2)/2,tot 



02 

'^2,(n-l)/2,tot 



02 

'-^2,-l/2,tot 



'2,(n-l)/2,tot'52,-l/2,tot 



;d.78) 



where I have used Eq.(D.74). 

Finally, if one considers the reconstruction only in bins, the diagonal entries of the 
covariance matrix in Eq.(4.49) can be reduced to 



[flA^s,n)) = ^ +VV 



F\Qs,n) 



H 2 



flri'Vs,n) 



(D.79) 



since r„ and kn are now two independent variables, and, similar to Eq.(5.31), I have 
defined here 

-1 



fi . n. — 



dQ 



In F^Q) 



Q=Qs 



(D.80) 



The expression in Eq.(D.79) is essentially the same as that in Eq.(D.72) without the 
last three terms involving <S'2,_i/2,tot) which correspond to the statistical error of the 
estimator for M and have been neglected in Eq.(4.49). Note that /i(w„) in Eq.(D.68) is 
obtained from the first expression of the exponential ansatz in Eq.(4.18) and estimated at 
V = v„, while /i .,,(fs,/i) in Eq.(4.46) is obtained from the second expression in Eq.(4.18) 
and estimated at f = Vs^^, which is not a fixed value like but actually depends on 
kn through Eqs.(4.28). Hence, since the two expressions in Eq.(4.18) are equivalent, 
when one takes into account the uncertainty of the determination of Qs,^ by Eq.(4.28), 
the statistical error of fir{vs,fj.) will be identical to the first two terms of cr^(/i(t>„)) in 
Eq.(D.72). 

Similarly, if one neglects Qthre and thus all terms involving rthre; the diagonal entries 
of the general form of the covariance matrix given in Eq.(4.61) can be reduced to 



,n\2„2. 



a\lo) + a^'^in + 1) V^(7„) - 2a"{n + l)(^;")cov(7o, 4) 



a\Io) ^ a\In) 2cov(J„,/o) 



-'0 



(D.81) 



where 1 have used (f") in Eq.(4.66). Comparing the definition of S'4^A,tot in Eq.(D.76) 
with the expression of cov(/„, in Eq.(4.54) and using Eq.(D.67), the expression given 
in Eq.(D.78) is exactly identical to that in Eq.(D.81). 
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Appendix E 

Some Detailed Calculations 



In this chapter I give some detailed derivations. 



E.l Derivations of covariances in Sec. 4.2 



E.1.1 Covariances in Subsec. 4.2.2 



Prom Eq.(4.37), one has 



1 



(E.l) 



where n* denotes a given bin between the n^_-th and the n^+-th bins. Then, from (4.38), 
one has 



N,.. 



and 



Combining Eqs.(E.l) to (E.3), it can be found that 
COY (Q - Q^,\f,,Q - QX) 



E 



dNn* 



+ 



dQ - Qi,\A [dQ-Q, 



dQ\ 



dQ\r 



a' (gin*) 



N N 



E (gin* - gu) (gin* - gi.) + N^a^ (g|, 



Meanwhile, from Eqs.(4.37) and (4.40), one has 



1 



W 



(E.2) 



(E.3) 



(4.39) 



(E.4) 



137 



Thus, 

dr^ ^ 1 
and then 



oov(v^J = E(|5-)(^)^=(iV..)^^i;iV„.. 



Combining Eqs.(E.2) and (E.5), it can be found that 

GOV (r^,g-g.|.) = 



dQ-Q, 



" ^2 



E.1.2 Covariance in Eq.(4.61) 

The expression of (v") in Eq.(4.52) can be rewritten as 

^Vthre ' thre 



i^2(Qthre) 



+ (n + 1)7„ 



with 



+ /o 



Hence, it can be found that 



A4.a"(n + 1), 



and 



Moreover, 



— —a 



• ' m 



^ Vthrc ' thrc 
i^'(Qthre) 



+ (n + 1)4 



dr 



thre 



^^(n+l)/2 
i^'(Qthre) 

2 



^Vthre ^thr 
i^'(Qthre) 



+ (n + 1)4 



(E.5) 



(4.41) 



(4.42) 



(4.52') 



(4.62') 



(E.6) 



(E.7) 



i^'(Qthre) 



(E.8) 



.i^'(Qthre). 

This leads to the definition of D„ in Eq.(4.63). Combining Eqs.(E.6), (E.7), and (4.63), 
Eq.(4.61) can be obtained directly. 
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E.2 Derivation of the correction terms in Eq.(4.64) 



Starting point is the observation that we wish to compute the ratio of two integrals, 

Gi ^ J gi{x)dx ^ J2i nigijxi) 
G2 Jg2{x)dx Ejnjg2{xj)' 

In the second step the integrals have been discretized, i.e., replaced by sums over bins i 
with rii events per bin. rii can be written as sum of average value rij and fluctuation dnf 

Gi J2i{ni + Sni)gi{xi) 



G2 J2j njg2{xj) + Ej Snjg2{xj) ' 
Introducing the notation 

Ga = Yl'^idaiXi) , 

i 

for a = 1, 2, and expanding up to second order in the 6ni, one has 
Gi Gi + Y.i5nigi{xi) 



(E.IO) 



(E.ll) 



G2 



G, 



^ _ Snjg2{xj) ^ I Ei Snjg2{xj] 



G, 



G, 



£1 + 
G2 G 



(j2Snigi{xi)^ - p (^Snig2{xi)^ 



(E.12) 



• 2 \ I / \ J y ^2 

Now consider the average over many experiments. Of course, 5ni averages to zero, but 
the product SriiSrij averages to fiiSij, i.e., it is non-zero for i — j. Hence: 



\G2 



£1 
G2 



=2 \J2^i9l{Xi)92{Xi)] + S (j2^i9l 

Go \ i I Gn \ i 



Xi 



(E.13) 



'2 \ ? J ^2 

The sums appearing in the two correction terms also appear in the definition of the 
covariance matrix between G\ and 6*2. Note that we wish to compute the first term on 
the right-hand side, since in this case the estimators for G\ and G2 indeed average to the 
correct values. This then leads to the final result 



Applying this result to Eq.(4.52) then immediately leads to Eq.(4.64). 

E.3 Statistical error of 7^n(<3thre) in Eq.(4.72) 



:e.i4) 



Prom Eq.(4.72), it can be found directly that 



5^n(Qthre) 



dr 



thre,X 



2 

n 



1/2 



X 



thrcX J 



^n(Qthre), (E.15a) 
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dTZniQthre) _n + l 



n 



thrcX J 



^^thTe,X 'thre,X 



in 



^)In,xFx{QthYe,x) 



^n(Qthre) j 



and 



dTZn{Qthie) 

dIo,x 



1 

n 



Fx{Qthre,x) 



T^n{Q thre) ■ 



(E.15b) 



(E.15c) 



_2Qjhre,x'^thre,X + Io,xF^(Qthre,x) _ 

By first exchanging Qthrc^x^ ^'^'^ (^+ ^)^n,x with Qthre.x /o,x, respectively, and then 
replacing X hj Y, one can get 



dTZniQthTe) 
drthre,Y 



2 



.2Qthi^,y'^^thre,y + ('T' + l)-^n,y-^y (Qthre,y) 

^y(Qthre,y) 

2Qthre,y'^thre,y + -^0,y-^y (Qthre,y) 



97?.„((5thre) 

9/n,y 



n + 1 



FyiQihrex) 



2QlCj n^re,Y + (n + 1) 

-^n,y-^y(Qthre,y) 



^n(Qthre), (E.16a) 



^n(Qthre) ) 

(E.16b) 



and 



dTZniQthie) 

dIo,Y 



1 



-^y(Qthre,y) 



2Qthre,y'^thre,y + -^0,y-^y (Qthre,y ) 



^n(Qthre) • 



(E.16c) 



By setting Qthre = 0, the above expressions can be reduced to cr(m^) given in Eq.(4.70) 
directly. 



E.4 Derivation of rjn in Eq.(5.24) 

Prom Eq.(5.22), I have 

/i,m(Q) ^ Q (dR\ j^_(dF\_ 1 



(m) 



M. 



(E.17) 



Substituting the ansatz for [dR/ dQ)(i)^n in Eq.(5.13), it can be found that 
1 



{dR/dQ)^l),n 

1 



(dR/dQ) 



(0) 



dQ V<^<5/(o) 



^ni,Q Qn) ~l~ 



(E.18) 



140 



Thus I can obtain that 

/l,(m)(Q) 



^ni,Q Qn) ~l~ 



X 



[ 2 /dF\ _ 1 
\F{Q)[dQ) (dR/dQ) 



(0) 



d /dR\ 
dQ [dQj 



(0). 



ln{Q - Qn) + K 



/l,(0)(Q) 

2J\f 



ln{Q - Qn) + K 



Q (dR\ 



(E.19) 



namely, 



/l,(l),n(Q) 

hmiQ) 

^n{Q Qn) ~l~ 



[ 2 



"|F(Q) [dQj {dR/dQ)^o) 
Now I use the ansatz for {dR/ dQ){iS)^n in Eq.(4.18) to get 



d (dR\ 

dQ VQl^^y 



-1 



. (E.20) 



/l,(l),n(Q) 
/l,(0),n(<3) 



ln{Q - Qn) + K 



F{Q) \dQ 



1 -1 



(E.21) 



Let Q = Qn: the expression of rjn in Eq.(5.24) can be obtained directly. 



141 



Bibliography 



[1] G. Jungman, M. Kamionkowski, and K. Griest, Phys. Rep. 267, 195 (1996). 
[2] F. Zwicky, Helv. Phys. Acta 6, 110 (1933). 
[3] S. Smith, Astrophys. J. 83, 23 (1936). 

[4] V. C. Rubin and W. K. Ford, Astrophys. J. 159, 379 (1970). 

[5] S. M. Faber and J. S. Gallagher, Annu. Rev. Astron. Astrophys. 17, 135 (1979). 

[6] V. C. Rubin, W. K. Ford, and N. Thonnard, Astrophys. J. 238, 471 (1980). 

[7] K. G. Begeman, A. H. Broeils, and R. H. Sanders, Mon. Not. R. Astron. Soc. 249, 
523 (1991). 

[8] R. R Oiling and M. R. Merrifield, Mon. Not. R. Astron. Soc. 311, 361 (2000). 
[9] M. Fich and S. Tremaine, Annu. Rev. Astron. Astrophys. 29, 409 (1991). 
[10] H. Stocker, Taschenbuch der Physik, Harri Deutsch (2005). 

[11] Review of Particle Physics, J. Phys. G: Nucl. Part. Phys. 33, 2. Astrophysical Con- 
stants and Parameters (2006). 

[12] O. Lahav and A. R. Lindle, Review of Particle Physics, J. Phys. G: Nucl. Part. 
Phys. 33, 21. The Cosmological Parameters (2006). 

[13] WMAP CoUab., D. N. Spergel et ai, arXiv:astro-ph/0603449 (2006). 

[14] E. W. Kolb and M. S. Turner, The Early Universe, Addison- Wesley (1990). 

[15] M. S. Turner, arXiv:astro-ph/9904051 (1999). 

[16] K. A. Ohve and J. A. Peacock, Review of Particle Physics, J. Phys. G: Nucl. Part. 
Phys. 33, 19. Big-Bang Cosmology (2006). 

[17] B. D. Fields and S. Sarkar, Review of Particle Physics, J. Phys. G: Nucl. Part. Phys. 
33, 20. Big-Bang Nucleosynthesis (2006). 



142 



[18] D. Scott and G. F. Smoot, Review of Particle Physics, J. Phys. G: Nucl. Part. Phys. 
33, 23. Cosmic Microwave Background (2006). 

[19] J. R. Primack, Nucl. Phys. Proc. Suppl. 124, 3 (2003). 

[20] J. A. Tyson, Phys. Scr. T 85, 259 (2000). 

[21] P. J. E. Peebles, Rev. Mod. Phys. 75, 559 (2003). 

[22] M. de Jesus, Int. J. Mod. Phys. A 19, 1142 (2004). 

[23] C. S. Frenk, Phil. Trans. Roy. Soc. Lond. 300, 1277 (2002). 

[24] M. Kamionkowski and A. Kinkhabwala, Phys. Rev. D 57, 3256 (1998). 

[25] Review of Particle Physics, J. Phys. G: Nucl. Part. Phys. 33, 1. Physical Constants 
(2006). 

[26] A. M. Green, Phys. Rev. D 63, 043005 (2001). 

[27] P. D. Sackett et ai, Astrophys. J. 436, 629 (1994). 

[28] N. W. Evans, Mon. Not. R. Astron. Soc. 267, 333 (1994). 

[29] J. J. Binney, Mon. Not. R. Astron. Soc. 196, 455 (1981). 

[30] N. W. Evans, Mon. Not. R. Astron. Soc. 260, 191 (1993). 

[31] J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 462, 563 (1996). 

[32] A. Burkert and J. Silk, arXiv:astro-ph/9904159 (1999). 

[33] G. Gentile, C. Tonini, and P. Salucci, arXiv: astro-ph/0701550 (2007). 

[34] E. L. Lokas, Mon. Not. R. Astron. Soc. 311, 423 (2000). 

[35] T. Kuwabara, A. Taruya, and Y. Suto, Publ. Astron. Soc. Jap. 54, 503 (2002). 
[36] L. Hernquist, Astrophys. J. 356, 359 (1990). 

[37] M. Drees and G. Gerbier, Review of Particle Physics, J. Phys. G: Nucl. Part. Phys. 
33, 22. Dark Matter (2006). 

[38] M. W. Goodman and E. Witten, Phys. Rev. D 31, 3059 (1985). 

[39] I. Wassermann, Phys. Rev. D 33, 2071 (1986). 

[40] K. Griest, Phys. Rev. D 38, 2357 (1988). 

[41] A. K. Drukier, K. Freese, and D. N. Spergel, Phys. Rev. D 33, 3495 (1986). 

143 



[42] K. Preese, J. Prieman, and A. Gould, Phys. Rev. D 37, 3388 (1988). 

[43] D. N. Spergel, Phys. Rev. D 37, 1353 (1988). 

[44] P. P. Smith and J. D. Lewin, Astropart. Phys. 6, 87 (1996). 

[45] Y. Ramachers, Nucl. Phys. Proc. Suppl. 118, 341 (2003). 

[46] CDMS CoUab., D. S. Akerib et al, arXiv:astro-ph/0609189 (2006). 

[47] EDELWEISS CoUab., V. Sanglard et al, Phys. Rev. D 71, 122002 (2005). 

[48] See http : //people . roma2 . inf n . it/~dama/web/home . html. 

[49] S. P. Ahlen et al, Phys. Lett. B 195, 603 (1987). 

[50] J. Engel, Phys. Lett. B 264, 114 (1991). 

[51] J. Gascon, arXiv: astro-ph/0504241 (2005). 

[52] DAMA CoUab., R. Bernabei et al, Phys. Lett. B 480, 23 (2000). 

[53] DAMA CoUab., R. Bernabei et a/., arXiv: astro-ph/0305542 (2003). 

[54] DAMA CoUab., R. Bernabei et al, arXiv: astro-ph/0311046 (2003). 

[55] See http : / /hepwww . rl . ac . uk/UKDMC/pro j ect/pro j act . html. 

[56] See http : //www . cresst . de/. 

[57] CRESST CoUab., M. Altmann et al, arXiv: astro-ph/0106314 (2001). 

[58] See http : //www . mpi-hd . mpg . de/ non_acc/dm . html. 

[59] See http : / /cdms . berkeley . edu/. 

[60] See http : //edelweiss . in2p3 . f r/index_newe . html. 

[61] L. Kaufmann and A. Rubbia, J. Phys. Conf. Ser. 60, 264 (2007). 

[62] R. W. Schnee, AIP Conf. Proc. 903, 8 (2007). 

[63] CDMS CoUab., D. S. Akerib et al, Phys. Rev. D 72, 052009 (2005). 
[64] CDMS CoUab., P. L. Brink et al, arXiv:astro-ph/0503583 (2005). 
[65] CDMS CoUab., D. S. Akerib et al, Phys. Rev. Lett. 96, 011302 (2006). 
[66] CDMS CoUab., D. S. Akerib et al, Phys. Rev. D 73, 011102 (2006). 

[67] CRESST CoUab., G. Angloher et al, Astropart. Phys. 23, 325 (2005). 

144 



[68] CDMS CoUab., R. Abusaidi et al, Phys. Rev. Lett 84, 5699 (2000). 

[69] P. Gondolo and G. Gelmini, Phys. Rev. D 71, 123520 (2005). 

[70] D. Tucker-Smith and N. Weiner, Phys. Rev. D 72, 063509 (2005). 

[71] V. Sanglard, for EDELWEISS CoUab., arXiv:astro-ph/0612207 (2006). 

[72] See http : //dmrc . snu . ac . kr/. 

[73] KIMS CoUab., H. S. Lee et a/., arXiv: 0704. 0423 [astro-ph] (2007). 

[74] PICO-LON CoUab., K. Fushimi et al, J. Phys. Soc. Jap. 74, 3117 (2005). 

[75] PICO-LON CoUab., K. Fushimi et a/., arXiv:astro-ph/0611700 (2006). 

[76] PICO-LON CoUab., K. Fushimi et al, J. Phys. Soc. Jap. 75, 064201 (2006). 

[77] WARP CoUab., P. Benetti et al, arXiv: astro-ph/0701286 (2007). 

[78] ZEPLIN CoUab., G. J. Alner et al, arXiv: astro-ph/0701858 (2007). 

[79] See http : //neutrino . ethz . ch/ArDM/. 

[80] A. Rubbia, J. Phys. Conf. Ser. 39, 129 (2006). 

[81] L. Kaufmann and A. Rubbia, arXiv:hep-ph/0612056 (2006). 

[82] See http://warp.lngs.infn.it/. 

[83] See http : / /xenon . astro . Columbia . edu/. 

[84] K. Ni and L. Baudis, arXiv:astro-ph/0611124 (2006). 

[85] XENON CoUab., J. Angle a/., arXiv: 0706. 0039 [astro-ph] (2007). 

[86] L. Baudis, arXiv:astro-ph/0703183 (2007). 

[87] N. Spooner, arXiv: 0705. 3345 [astro-ph] (2007). 

[88] See http : //Ipse . in2p3 . f r/inimache3/. 

[89] E. Moulin and D. Santos, arXiv:astro-ph/0505458 (2005). 

[90] D. Santos et al, arXiv:astro-ph/0701230 (2007). 

[91] SIMPLE CoUab., T. Girard et al, arXiv:hep-ex/0504022 (2005). 

[92] SIMPLE CoUab., T. Morlet aZ., arXiv: 0704. 2037 [astro-ph] (2007). 

[93] SIMPLE CoUab., T. Girard et al, Phys. Lett. B 621, 233 (2005). 

145 



[94] Review of Particle Physics, J. Phys. G: Nucl. Pari. Phys. 33, 3. International System 
of Units (SI) (2006). 

[95] F. S. Ling, P. Sikivie, and S. Wick, Phys. Rev. D 70, 123503 (2004). 

[96] M. Drees and C. L. Shan, J. Cosmol. AstropaH. Phys. 0706, Oil (2007), 
arXiv : astro-ph/0703651. 

[97] A. M. Green, arXiv:hep-ph/0703217 (2007). 



146 



